{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Modeling and Simulation in Python\n",
    "\n",
    "Insulin minimal model\n",
    "\n",
    "Copyright 2017 Allen Downey\n",
    "\n",
    "License: [Creative Commons Attribution 4.0 International](https://creativecommons.org/licenses/by/4.0)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Configure Jupyter so figures appear in the notebook\n",
    "%matplotlib inline\n",
    "\n",
    "# Configure Jupyter to display the assigned value after an assignment\n",
    "%config InteractiveShell.ast_node_interactivity='last_expr_or_assign'\n",
    "\n",
    "# import functions from the modsim.py module\n",
    "from modsim import *"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Data\n",
    "\n",
    "We have data from Pacini and Bergman (1986), \"MINMOD: a computer program to calculate insulin sensitivity and pancreatic responsivity from the frequently sampled intravenous glucose tolerance test\", *Computer Methods and Programs in Biomedicine*, 23: 113-122.."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "data = pd.read_csv('data/glucose_insulin.csv', index_col='time');"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### The insulin minimal model\n",
    "\n",
    "In addition to the glucose minimal mode, Pacini and Bergman present an insulin minimal model, in which the concentration of insulin, $I$, is governed by this differential equation:\n",
    "\n",
    "$ \\frac{dI}{dt} = -k I(t) + \\gamma (G(t) - G_T) t $"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Exercise:**  Write a version of `make_system` that takes the parameters of this model, `I0`, `k`, `gamma`, and `G_T` as parameters, along with a `DataFrame` containing the measurements, and returns a `System` object suitable for use with `run_simulation` or `run_odeint`.\n",
    "\n",
    "Use it to make a `System` object with the following parameters:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>values</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>I0</th>\n",
       "      <td>360.000</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>k</th>\n",
       "      <td>0.250</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>gamma</th>\n",
       "      <td>0.004</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>G_T</th>\n",
       "      <td>80.000</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "I0       360.000\n",
       "k          0.250\n",
       "gamma      0.004\n",
       "G_T       80.000\n",
       "dtype: float64"
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "params = Params(I0 = 360,\n",
    "                k = 0.25,\n",
    "                gamma = 0.004,\n",
    "                G_T = 80)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Solution\n",
    "\n",
    "def make_system(params, data):\n",
    "    # params might be a Params object or an array,\n",
    "    # so we have to unpack it like this\n",
    "    I0, k, gamma, G_T = params\n",
    "    \n",
    "    init = State(I=I0)\n",
    "    \n",
    "    t_0 = get_first_label(data)\n",
    "    t_end = get_last_label(data)\n",
    "    G=interpolate(data.glucose)\n",
    "    \n",
    "    system = System(I0=I0, k=k, gamma=gamma, G_T=G_T, G=G,\n",
    "                    init=init, t_0=t_0, t_end=t_end, dt=1)\n",
    "\n",
    "    return system"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>values</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>I0</th>\n",
       "      <td>360</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>k</th>\n",
       "      <td>0.25</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>gamma</th>\n",
       "      <td>0.004</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>G_T</th>\n",
       "      <td>80</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>G</th>\n",
       "      <td>&lt;function interpolate.&lt;locals&gt;.wrapper at 0x7f...</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>init</th>\n",
       "      <td>I    360.0\n",
       "dtype: float64</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>t_0</th>\n",
       "      <td>0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>t_end</th>\n",
       "      <td>182</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>dt</th>\n",
       "      <td>1</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "I0                                                     360\n",
       "k                                                     0.25\n",
       "gamma                                                0.004\n",
       "G_T                                                     80\n",
       "G        <function interpolate.<locals>.wrapper at 0x7f...\n",
       "init                             I    360.0\n",
       "dtype: float64\n",
       "t_0                                                      0\n",
       "t_end                                                  182\n",
       "dt                                                       1\n",
       "dtype: object"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Solution\n",
    "\n",
    "system = make_system(params, data)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Exercise:** Write a slope function that takes state, t, system as parameters and returns the derivative of `I` with respect to time.  Test your function with the initial condition $I(0)=360$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Solution\n",
    "\n",
    "def slope_func(state, t, system):\n",
    "    [I] = state\n",
    "    k, gamma = system.k, system.gamma\n",
    "    G, G_T = system.G, system.G_T\n",
    "        \n",
    "    dIdt = -k * I + gamma * (G(t) - G_T) * t\n",
    "    \n",
    "    return [dIdt]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[-90.0]"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Solution\n",
    "\n",
    "slope_func(system.init, system.t_0, system)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Exercise:** Run `run_ode_solver` with your `System` object and slope function, and plot the results, along with the measured insulin levels."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>values</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>success</th>\n",
       "      <td>True</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>message</th>\n",
       "      <td>The solver successfully reached the end of the...</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "success                                                 True\n",
       "message    The solver successfully reached the end of the...\n",
       "dtype: object"
      ]
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Solution\n",
    "\n",
    "results, details = run_ode_solver(system, slope_func)\n",
    "details"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>I</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>178</th>\n",
       "      <td>22.356671</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>179</th>\n",
       "      <td>23.180483</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>180</th>\n",
       "      <td>24.013211</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>181</th>\n",
       "      <td>24.854654</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>182</th>\n",
       "      <td>25.704657</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "             I\n",
       "178  22.356671\n",
       "179  23.180483\n",
       "180  24.013211\n",
       "181  24.854654\n",
       "182  25.704657"
      ]
     },
     "execution_count": 9,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Solution\n",
    "\n",
    "results.tail()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAagAAAEYCAYAAAAJeGK1AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjAsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+17YcXAAAgAElEQVR4nO3deXzU9bX4/9fMJJOdBAIhgIhI4UjQWBC5LnWrdnG31drrQrW299qq1ZbaKtWr1VbRVimK9kqr1eoP9fYrFRTbutUFalVwQwQPArIvAbKSfZnfH5/PhMmQwMxktpDzfDzyyMxnPXwy5OT9/rw/7+MJBAIYY4wx6cab6gCMMcaY7liCMsYYk5YsQRljjElLlqCMMcakpYxUB5BsIpIFHA1sBdpTHI4xxvR3PmAYsERVm0NX9LsEhZOcFqU6CGOMMV2cACwOXdAfE9RWgLlz51JaWprqWIwxpl/btm0bl1xyCbi/m0P1xwTVDlBaWspBBx2U6liMMcY49rrlYoMkjDHGpCVLUMYYY9KSJShjjDFpyRKUMcaYtNQfB0kYY/qQ2tpaKioqaG1tTXUoJgaZmZmUlJQwYMCAqPe1BGWMSVu1tbVs376dESNGkJOTg8fjSXVIJgqBQIDGxkY2b94MEHWSsi6+KFU31bKzoTLVYRjTL1RUVDBixAhyc3MtOfVBHo+H3NxcRowYQUVFRdT7W4KK0q2v3ssNL95JW4fNkmRMorW2tpKTk5PqMEwv5eTkxNRFawkqSrUtu6lrqaextTHVoRjTL1jLqe+L9WdoCSpKWT4/AM3tLSmOxBhjIrNx48ZUhxCTlCYoETlLRJaJSJ2IrBGRK93lWSLSIiK7Q75eCtnvQhFZLSL1IvI3ESlJVszBBNXSZgnKGLN/Dz30ENOmTYvrMTdt2oSIUFtbu99t586dy1133dX5fuLEiahqXONJlJSN4hORYcAzwDdU9e8iMgn4l4gswUmclaq612yuIlIGPAKcDiwF7gaeBr6cjLizMpwE1WQJyhgTgR/84AcpPX9lZSWBQKDz/QcffJDCaKKTshaUqm4FhrjJyQsUA21AHXAU8GEPu14KPK+qi1W1CZgOHC8iY5MRd2cLyrr4jDEhOjo6uPPOOznuuOM49thj+d73vseGDRuYPXs2V111FQCzZ89m+vTpXHXVVUycOJGzzz6bDz/8kGuvvbbzfbB1E7of7LvV9OKLL3L++eczZcoUjj76aKZPn05raysvvvgic+bM4fXXX+ecc84BQERYuXIlAMuXL2fq1KlMnjyZr33ta8ydO7fzmFOnTuV3v/sd3/jGN5g0aRIXX3wxa9asSdj1605Kn4NS1ToRyQVq3FjuVtXPROR6oERElgFDgTeBH6vqZqAMp+UUPEaDiGwEjgA+S3TM/gy7B2VMKs1480E+2Lo8KeeaOOxwpp94dUTbvvzyy7z55pv8/e9/Jzc3l1tuuYU5c+bsVdbnueee48EHH2T27NlcffXVXHrppTz44IPcc8893HjjjTzwwAPMnj074hg3b97Mz3/+cx555BEmT57M+vXr+fa3v80rr7zC6aefzqpVq1i5ciW///3vu+xXWVnJ5ZdfztVXX82f/vQnVq1axZVXXklhYSFnnXUWAPPnz+exxx5jyJAhXHfddcyePZtZs2ZFHFtvpcMgiSYgD6eQ4BUi8j2gHvgXcCogQCPwrLt9PtAQdowGIDcZwXYOkrAuPmNMiIKCAioqKliwYAHbtm3jjjvu4I477thru/Lyck4++WR8Ph9TpkxhzJgxnHTSSfj9fo477jg2bdoU1XmHDBnCwoULmTx5MnV1dVRWVjJw4MD9Pnf06quvMmTIEL773e+SmZnJhAkT+M53vsO8efM6tznnnHMYPXo0+fn5fO1rX2PDhg1RxdZbKZ9JQlU7gBZgqYj8AThXVc8J3UZEpgE7RGQkTvIKfzAiF9idjHiDLSjr4jMmNSJt0STbcccdxy233MLTTz/NXXfdxciRI7nhhhv22q6oqKjztdfr7TK7gtfrpaOjI6rzZmZmMm/ePJ555hmys7MpKyujubm5y32n7lRWVjJ8+PAuy0aMGMHWrXvqBhYXF3e+zsjIoL09uc9/pqwFJSInich7YYuzgGoRuV1Exocs97vfm4AVOK2q4HFygYPd5QlnLShjTHc2btxIWVkZTz/9NO+88w7f/OY3+fGPf0xbW1uX7SJ9Jsjr9XZ5uLW6urrb7V544QWef/555s2bxyuvvML9999Pfn7+fo8/bNgwtmzZste/YfDgwRHFlwyp7OL7EBghItNExCcixwHfwxmhVw7cKyJFIlIE3Ae8oKo7gCeBc0XkZBHJAmYAH6jqqmQEnWX3oIwx3Xj77be55ppr2LJlC/n5+RQWFlJQUIDP54vpeKNHj+aDDz5g7dq1NDQ08Oijj3a7XV1dHT6fD7/fT2trK0888QSq2pnc/H4/dXV1e+130kknUVVVxWOPPUZraysrVqzgiSee4Oyzz44p3kRI5Si+GuAM4JtAJfAH4Puq+gZOoqoCVgPrcLoAp7r7fQxcATwE7AQmAN9KVtzWgjLGdOf888/ntNNO48ILL2TSpEn85S9/4f777495FoXTTjuNM888k4suuogzzjiDY489ttvtvvGNb1BWVsZpp53GiSeeyNtvv81ZZ53FZ585Y8ZOPvlk1q9fz0knndRlv8LCQh5++GFeeeUVjjnmGK655hq+//3v8+1vfzumeBPBs79+ygONiBwCfP7qq69y0EEHRb3/M5+8wF+WL+SbZafzn0ecs/8djDExW7lyJePHj9//hibt9fSz3LRpE6eeeirAaFVdF7ouHUbx9Sl+m0nCGGOSwhJUlGwuPmOMSQ5LUFHqHCRhLShjjEkoS1BRslF8xhiTHJagomRz8RljTHJYgoqS34aZG2NMUliCipJ18RljTHJYgoqSFSw0xpjksAQVpc6ChdaCMsb0I+3t7V0mkk0GS1BRshaUMSbcli1bmDhxYrdz3sXL1KlTeeyxxwD4/ve/36W4YKz++te/cu6550a07bRp03jxxRd7fc5opLzcRl9jBQuN6VsWrX+Xp5YtYFdDJcW5g7io/FxOGDUlrucYPnx4UkupP/zww0k7V1BlZWXSz2ktqCj5Q2aS6G/zGBrT1yxa/y5zlsxlZ0MlAWBnQyVzlsxl0fp343qe0HLsmzZtYuLEiTz66KN86Utf4thjj+XWW2/trPP0+uuvc+aZZzJ58mTOPvts5s+fv9cxgkJbTaFCl0dTmr26upof/ehHTJo0ia9//essX76nMnEgEOCBBx7g9NNPZ+LEiZx44omd57jjjjtYunQp99xzD7fffjsATz75JGeffTZHHXUUxx57LL/97W97exn3YgkqShleHz6vj0AgQFtH2/53MMakzFPLFuz1zGJLewtPLVuQ0PM2NDSgqrzyyis88sgjPPfccyxatIiOjg6uv/56fv7zn7N06VKmT5/Or3/9a+rr63t1vvnz5zNz5kzefPNNcnJyeiwZf8stt9DS0sKbb77JnDlzeOONNzrXLVy4kAULFvDYY4/x/vvvc+utt/Kb3/yGiooKbrrpJiZPnsz111/PLbfcwvvvv8+sWbOYNWsW7733HnPmzOHPf/4zy5Yt69W/I5wlqBjYfHzG9A27GrrvluppeTz993//d2eFWxFhw4YNeL1e8vLyeOGFF1i6dClTpkzh3XffJS8vr1fniqQ0e3NzM//85z/50Y9+RH5+PqNGjWLq1Kmd60855RTmzp3L0KFD2blzJ5mZmbS3t3fbtTd+/Hjmz5/PmDFjqKqqoqmpiby8vP2WmY+W3YOKQVaGn4bWRprbWsj39+6DZYxJnOLcQezsJhkV5w5K+LlDK9MGf9kD/OlPf2L27NlcddVVtLe3c/755/Ozn/2sV+eKpDR7dXU1ra2tlJaWdi4LLTnU1tbGjBkzeOuttygpKaG8vByg21sZPp+POXPm8OKLLzJw4EDKysqiLlUfCUtQMbAWlDF9w0Xl5zJnydwu3Xx+n5+LyiMbuRZvjY2NVFRUMGvWLDo6Onj//fe59tprmTBhAlOmOAM3IinzHouBAwfi9/vZsmVLZ/Lcvn175/qZM2fS3NzMG2+8QXZ2NjU1NTzzzDPdHuvRRx9lxYoVvPTSSwwYMIBAIMDRRx8dt1iDrIsvBjbU3Ji+4YRRU7jy6EsYnDsIDzA4dxBXHn1J3EfxRaq9vZ0f/vCHvPDCC3g8HkpLS/F4PBQVFVFcXExBQQHz58+nvb2dv//97z0OdoiF3+/nrLPOYtasWdTU1LBp0yYef/zxzvV1dXVkZWXh8/moqanhzjvvBOhSOn737t2d22ZmZpKRkUFjYyMzZ86krq6Olpb4/k60FlQMbKi5MX3HCaOmpCwhhcvPz+f+++/nnnvu4eabbyY/P59LLrmksxz7nXfeyb333suDDz7IKaecwle+8pW4nv/mm2/ml7/8JaeccgpFRUWcdtppvPPOOwBcd9113HjjjUyZMoWCggLOOOMMRIRVq1ZRXl7O2Wefze23387nn3/OTTfdxMqVKzn++OPJzc3lxBNP5Pjjj+8sMx8vKS35LiJnAXcCo4EK4DeqOkdE/MADwAVAOzBTVWeE7Hehu98w4A3gclWN6O5cb0u+A9z+2iyWVyg3n3Qt5aVWjtqYRLGS7weOPlXyXUSGAc8AN6hqAfAtYJaITAJuAwQYAxwNXCYi33H3KwMeAS4HioHPgKeTGbu1oIwxJvGi7uITES9wGFCC07rZBqxW1aiaYqq6VUSGqGqde8xioA2oAy7DaRVVAVUicg9wJfA4cCnwvKouduOZ7m4zVlXj277sQbaV3DDGmISLOEGJyInAdcBpQEHIqgBOgvgH8HtVfSvSY7rJKReocWO5G9iB03W3ImTTT4Ej3NdlwNKQYzSIyEZ3fVISVLAFZUULjTEmcfaboERkLDAHOBh4FvgmTvLYhdNFOAQ4EjgReFpE1gBXquqqCGNoAvKAcuBvQKO7vCFkmwYg132dH7YufH3CZVkLyhhjEi6SFtT/B9yuqi/0sH6j+7VQRG4AznP3iWjYjKp2AC3AUhH5AzDZXZUTslkusNt9XR+2Lnx9wlnRQmOSJxAI4PF4Uh2G6YVYB+NFMkjimH0kpy5UNaCqzwL/sb9tReQkEXkvbHEWUIVzX0tClh/Gni6/FaHr3C7Cg+naJZhQwQljrYvPmMTKzMyksbFx/xuatNbY2EhmZmbU++23BRXJ4AcRGQ6cpqqPR7oP8CEwQkSmAffhJLXvAd/ASVC3isgynC69691tAJ4EFovIycC/gRnAB1F0KfZadrBooXXxGZNQJSUlbN68mREjRpCTk2MtqT4mEAjQ2NjI5s2bGTp0aNT7x+tB3SOAR3FG2UVEVWtE5AzgfuBWnG7C76vqGyLyDnAv8AlOK+8PwEPufh+LyBXu+xHAOzhD1JPGbzNJGJMUAwYMAJyCgKFTAJm+IzMzk6FDh3b+LKOR0pkkVPV94EvdLG8Crna/uttvHjAvsdH1zObiMyZ5BgwYENMvN9P32Vx8MbBBEsYYk3iWoGJgXXzGGJN4kTwHdUYExzkqDrH0GdaCMsaYxIvkHtTCCI+Vullnk8we1DXGmMSLZJi5dQOG6WxBtTWnOBJjjDlwxTJZbC5OeYyssFUBVf0gLlGluZyMbAAaLUEZY0zCRJWgRORSnOePcoDwJ+YCgC9OcaW17MxggmpKcSTGGHPgirYFNQOnFtNMnEle+6WcDKfx2NTaZPOEGWNMgkSboAYAD6jq+kQE01dk+jLxeX20d7TT2tGG3xf9HFPGGGP2LdoBEE/gVLLt93KD96FabSJLY4xJhGhbUL8F3heRS4B1QEfoSlX9cpziSnvZmdnUtdTT2NZMYaqDMcaYA1C0CeoJnLpLL7B30cB+JTiSr6m1396KM8aYhIo2QR0N/IeqLktEMH1JcKCEjeQzxpjEiPYelAJFiQikr8kJDjW3FpQxxiRELMPMHxORB4A1QJcCLar6t3gFlu7sWShjjEmsaBPUU+73e7pZ128e1IWQ2SRabTYJY4xJhKgSlM3Lt0ewi6/JWlDGGJMQkZTb+Ax4CXgZeFVV6xIeVR8QbEE12D0oY4xJiEhaUGcCXwG+C/xJRFayJ2G9raod+9r5QJWTuWe6I2OMMfEXSbmNVcAq4EER8QHH4iSse4GxIrIIN2Gp6upoTi4iXwHuAsYCFcBvVXWOiGQBdUBowaW3VPWr7n4XAncCw4A3gMtVtSKac/dWts1obowxCRXtPah2YLH7dauIFAKn4iSsaTiJJiIiMhKYB1wGLMCpyvuiiKwDdgGVqlrazX5lOBPWng4sBe4GngaSOotFro3iM8aYhIrkHlRZD6tagWpV/Svw1xjOfQjwpKo+675fIiKvA8cDm4EPe9jvUuB5VV3sxjcdqBKRsar6WQxxxKSzBWVdfMYYkxCRtKCW4wwh766mREBEtgJ3qeoD0ZxYVRcBi4LvRWQQcALOdEpfB0pEZBkwFHgT+LGqbgbKcFpOweM0iMhG4AggaQnKRvEZY0xiRTJsfDRwqPs99GsMztRHdwO3icjlsQbhdhU+B7yD091XD/wLp/tQgEYg2NLKZ+95ABuA3FjPH4vOqY6sBWWMMQkRySCJfdV++hxndvNG4DrgsWgDEJFxOElpBXCJOypwWtg204Ad7n2repyKvqFycSaxTZqcTCcEGyRhjDGJEa8Hb98kigESQSJyIk6raT5wgao2uctvF5HxIZv63e9NOIlMQo6RCxzsLk+aPS0oqwdljDGJEO1URz3xAlE1JURkDLAQuElVZ4etLgcmi8jF7vv7gBdUdYeIPAksFpGTgX/jzA/4gTscPmn2zMVnLShjjEmEeLWgLgfei3Kfq4ECYIaI7A75uhv4HlAFrMYpjNgCTAVQ1Y+BK4CHgJ3ABOBbcfg3RCXL58fj8dDa3kp7R3uyT2+MMQe8SIaZX9XDKi9QCByHM5jh1GhOrKrTCLvXFOaSfew7D+cZqpTxeDzkZGTT0NpIY1sT+f68VIZjjDEHnEi6+H7Ww/JWnFbOB8DRbsumX+lMUK2WoIwxJt4iGcU3OhmB9EXZmVnQaEPNjTEmEaIaJCEiPwDeBz5W1caQ5TkAocv6g9yM4MO6NlDCGGPiLdpRfNOBkUC7iChO9977QBtwLTEMNe/LrKquMcYkTlSj+FR1FDAYOANnSqIc4JfArGiPdSDIsfn4jDEmYaJ+DkpVK3FqQb0MICJFOJPF/m98Q0t/2Zk23ZExxiRKr1s9qloN/A9OfaZ+JTcjON2RJShjjIm3qBKUiBzUw6odQE/rDljBFpQNkjDGmPiLtotvg4jsxJk14n33azPOrBDPxDm2tBe8B9VgXXzGGBN30SaocmCi+3UccBXObBIA/xKR+3EKDX6oqu/HLco01VkTyhKUMcbEXbQl35fjFDB8IrhMREazJ2lNBM4HSgFf/MJMT7luyY361vDyVMYYY3orkrn4ityBEN1S1c9x6kL9NWSfkviEl97y/U6NxPoWS1DGGBNvkQySeENEbnCr3u6TiAwWkZuAV3sfWvrLswRljDEJE0kX3/HAr4FNIvIv4B/AJzilLjzAEOBI4CTgBODP7j4HvOAEsbstQRljTNxFMlnsbuDHInIXcCVwMc69puA9placKY9eAP5LVbckKNa0E2xB7bZ7UMYYE3cRD5JQ1W3AbcBtIuIFioEOVd2VqODSXX7mni6+QCCAx+NJcUTGGHPgiKnku6p24Dyc269l+DLIysiiua2ZxramzlF9xhhjeq/fTfAab6GtKGOMMfETUwsqXkTkK8BdOGU6KoDfquocEfEDDwAXAO3ATFWdEbLfhThz/w0D3gAuV9WKZMcPzn2oXY1V1Lc0MCSvOBUhGGPMASllLSgRGQnMwxkhWARcBMwQka/h3OsSYAxwNHCZiHzH3a8MeAS4HOc+2GfA08mOPyj4LNTulvpUhWCMMQekmFpQIpLp7ttlVICqRtPPdQjwpKo+675fIiKv4wxRvwynVVQFVInIPTgjCB8HLgWeV9XFbizT3W3Gqupnsfx7eqNzJJ918RljTFxFW/L9GGAOcHgPm0Q8vZGqLgIWhRx7EM5zVE/gdN2tCNn8U+AI93UZsDTkOA0istFdn7IEZfegjDEmvqJtQc0CaoDzgNp4BeHOUvEc8A7OTOkAob/xG4Bc93V+2Lrw9UllD+saY0xiRJugjgCOUdWP4xWAiIwDFuC0mC7BKSNPyHdwks9u93V92Lrw9UnVOR+fPaxrjDFxFe0giZXA8HidXEROxGk1zQcuUNUm977TNpxBEkGHsafLb0XoOhHJBQ6ma5dg0uS5w8x3N9sgCWOMiadoW1CzgT+KyGyc+z0toStV9W+RHkhExgALgZtUdXbY6ieAW0VkGU6X3vXAfe66J4HFInIy8G9gBvCBqq6K8t8SF/lZNt2RMcYkQrQJ6lH3+93drAsQXQ2oq4ECnKHlM0KWPwjcAtyLMymtF/gD8BCAqn4sIle470fgtMC+FcV54yov07kHZYMkjDEmvqItWBi356ZUdRowbR+bXO1+dbfvPJxnqFLOnoMyxpjEiPU5qFOBCTitm5XAq6raFs/A+gorWmiMMYkR7XNQpTgDGiYB63Ae1B0FfCoip6VquqFUsuegjDEmMaLtsrsPaANGq+o4VR2LMyNEJTAzzrH1CcFRfPWtjXR0dKQ4GmOMOXBEm6C+DlyrqpuDC9wChT8FzohnYH2F1+vtLLPR0NqY4miMMebAEW2CasIZrRcu2hF8BxSrrGuMMfEXbYJ6CZgpIkODC9zX9wIvxjOwviTfHtY1xpi4i3YU38+AfwLrRWS9u2wUsAy4OJ6B9SXBh3VtuiNjjImfaJ+D2iYi5cDXcIaZNwIrVfWVRATXV9jDusYYE3/7TVAicgbwsqq2uq/Buee03H3tDy6PZqqjA0l+lpOgaptTMl+tMcYckCJpQS0ESnFKsi/cx3b9dqBEYVYBALXNdSmOxBhjDhz7TVCh0xvFc6qjA0lhtpOgaposQRljTLxElXBE5J8iUtTN8iEi8l53+/QHlqCMMSb+IrkHdTJOmXWAk4ArRST8N/F4YEx8Q+s7CrMGAFBjXXzGGBM3kdyD2oVTj8njfl0NtIesD+BUs/1p3KPrI/a0oGpTHIkxxhw4IrkH9TFwKICIvAZ80616a1ydCcpaUMYYEzfRPgd1SnfLRcQPHKWq/45LVH1MXmYuPq+PxtYmWtpb8fsyUx2SMcb0edGW2zgGmINzTyp8gEUg2uMdKDweD4VZBVQ2VlPbVMfgvEGpDskYY/q8WMpt7AD+E2cWiUuB6Tj3oC6Kb2h9S/BZqGq7D2WMMXERbYIqB6a5Jdc/ALar6m+A64Afxzu4viR4H8oe1jXGmPiItkuuDQg2EVYBX8SZPPY1YFasQYjIFGChqpa477OAOqAlZLO3VPWr7voLgTuBYcAbwOWpruY7IDvYgrIEZYwx8RBtgnoHuEpEbgQ+As7GqaR7OF2TSURExAN8D7gnbNURQKWqlnazTxnwCHA6sBS4G3ga+HK054+nomznWShrQRljTHxE28U3HfguMA14AigTkbXAX4AnYzj/bcAPgV+HLT8K+LCHfS4FnlfVxara5MZ0vIiMjeH8cTPA7kEZY0xcRZWgVHUJcAjwZ/dZqKNwWj9XAD+J4fwPqepROC2hUJOAEhFZJiLbReT/icgId10ZsCIkpgZgI06rK2U6W1DWxWeMMXER7TDzhcDPVHUlOPWhgN/HenJV3dLDqnrgX8DtQCtwP/AsMAXIB8ILLzUAubHGEQ/BFlRNs7WgjDEmHqK9B3UMTsJIKFWdFvpeRKYBO0RkJE7yygnbJRdnqHvKFHVOd2Q1oYwxJh6iTVC/Ax4Xkd8Ba3Geheqkqiu63StKInI78FSwpQb43e9NON17ErJtLnAwId1+qVDodvHZfHzGGBMf0SaoX7nf/y9kWQBnEtl4FiwsByaLyMXu+/uAF1R1h4g8CSx2Z1n/NzAD+EBVV8Xp3DEpyMoHoLZlNx0dHXi9VjrLGGN6I9rfoqO7+To05Hu8fA+oAlYD63CGsE+FzslrrwAeAnYCE4BvxfHcMcnw+ijw5xEIBKhrsW4+Y4zprWhbUI/izGZeHbpQRIYA/8AZ1Rc1VX0dKAp5vwu4ZB/bzwPmxXKuRCrMHkBdSz3VTbWdXX7GGGNiE0vBwv8WkfAmQr8uWBhUnFvEptqt7GqoZlTRQakOxxhj+rRYChZegxUs7FZxzkAAdjZUpjgSY4zp+6xgYZwsWv8ub2/6AIC5Hz1LTmY2J4yakuKojDGm74qpYKGIZLr7esLWhz9A2y8sWv8uc5bMpaXdmY6wsa2JOUvmAliSMsaYGEU1ik9E/kNEPsJ5Hmk3zozjdSGv+6Wnli3oTE5BLe0tPLVsQYoiMsaYvi/aUXz3ATXAeewpu9Hv7erhnlNPy40xxuxftAnqCOAY976UcRXnDup2YERxrpV+N8aYWEX7oO5KYHgiAunLLio/F7/P32WZ35fJReXnpigiY4zp+6JtQc0G/igis4HPCCtSqKp/i1dgfUlwIMRTyxZ0tqTOG/91GyBhjDG9EMtMEuBUsQ0Xz7n4+pwTRk3hhFFTuHvR73lvy8eMLByW6pCMMaZPi3aYuc2Auh/Fuc7Dursa7FExY4zpjWhbUACISBlOyYuXgBJgnaoG4hlYXzXYHRix0xKUMcb0SrTPQQ0Qkb8By4FngKE4NaI+EhGbfA6b7sgYY+Il2i67e4Es4CD2lF2/FqgGZsUxrj5rcJ518RljTDxEm6DOBH6mqluCC1R1A/Aj4MvxDKyvCj77ZAnKGGN6J9oElU9YmXeXL4ZjHZAG5RTh9Xipaqyhpa1l/zsYY4zpVrRJ5R/AL93JYgECbrHCe4CX4xpZH5Xh9TE0bzABAmzbvSPV4RhjTJ8VbYL6ETASp0ZULvAKsBEoBK6Lb2h917CCEgC27q5IcSTGGNN3Rfsc1HbgOLfK7gR3/xXAK70ZZi4iU4CFqlrivvcDDwAX4BRHnKmqM0K2vxC4ExgGvAFcrqppkw2GFQyFrcvZWpc2ITjE9T0AABzuSURBVBljTJ8T7TDzDBG5HRBVfVBV7wN+DdwsIlHfgxIRj4h8H+d5qtDJ7G7Dec5qDHA0cJmIfMfdpwx4BLgcKMaZcunpaM+dSMEW1Ja67SmOxBhj+q5ok8o9wFRgbciyB3GSxa9iOP9twA9xklyoy4A7VLVKVde5573SXXcp8LyqLlbVJmA6cLyIjI3h/Akx3E1Q26wFZYwxMYs2QX0buFhVOwdEqOrjOAnl8hjO/5CqHgUsDS4QkSKcrrsVIdt9ilPqA6AsdJ1bxXdjyPqUG1YwFMC6+IwxpheiTVC5ONVzw1XhDJSISujzVCHy3e+h5eMb3HMH14eXlg9dn3IDcwrJ8vmpaa6jviU8VGOMMZGINkG9BvxGRDor8bktnjuA1+MUU737PSdkWWhirA9bF74+5bweL6XBkXzWijLGmJhEm6CuxRm4sFlEPhWRlcAW4FCcIei9pqpVwDacQRJBh7GnW29F6DoRyQUOpmuXYMoNswRljDG9Eu0w8w0icjjwFWA8TsHCVcBLqtoRx7ieAG4VkWU4XXrXA/e5654EFrtD3f8NzAA+UNVVcTx/rw3vfBbKRvIZY0wsoh4arqotOMPCHwQeBt4Est2WTLzcgjNj+ifAEmAe8JB7/o+BK9z3O3Gex/pWHM8dF8PynYESW2otQRljTCyiakGJyDHAHODwsFUeelFRV1VfB4pC3jcBV7tf3W0/Dydppa0RA0oB2FjT3TgQY4wx+xNtwcJZQA1wHlAb/3AOHBtrt7jft/LD537BxUeexwmjpqQ4KmOM6TuiTVBHAMe43WymB4vWv8sj7/1f5/tdjVXMWTIXwJKUMcZEKNp7UCuB4YkI5EDy1LIFtLR3LbXR0t7CU8sWpCgiY4zpe6JtQc0G/igis3HmwOvyW1hV/xavwPqyXT2Ue+9puTHGmL1Fm6Aedb/f3c26mAdJHGiKcwexs5tkFKy2a4wxZv+ifQ7KquZG4KLyc5mzZG6Xbj6/L5OLys9NYVTGGNO3RNuCAkBETsV5/siLc1/qVVVti2dgfVlwIMRTyxZ0tqS+Mf50GyBhjDFRiPY5qFJgPjAJWIfz/NMo4FMROS2digam2gmjpnDCqCnM+vcjvLVhKQNzBqQ6JGOM6VOi7bK7D2gDRqvqOFUdCxwCVAIz4xzbAWHMwFEA6M61+9nSGGNMqGgT1NeBa1V1c3CBWzLjp8AZ8QzsQDGhxKmj+EmFpjgSY4zpW6JNUE04o/XC2Qi+HhxSNJK8zBwq6ndRsXtnqsMxxpg+I9oE9RIwU0SGBhe4r+8FXoxnYAcKr9dLWck4AJZXpNWE68YYk9aiTVA/A4YC60VERUSB9UAecF28gztQHF7ilK9avv3TFEdijDF9R7TPQW0TkXKce1FlQCOwUlVfSURwB4rDh7oJqkIJBAJ4PJ4UR2SMMekvohaUiGSJyA9EZKCqtqnqQlX9DZDlrBZ/YsPs2w4aMIzC7AFUN9WyuXZbqsMxxpg+Yb8JSkQKcYoS/o6uZdgBhuFMe/RPESmIf3gHBo/HQ/nQwwBYumVZiqMxxpi+IZIW1P8AucBYVX07dIWq/hSneOEQ4BfxD+/AcezISQD8a8PSFEdijDF9QyQJ6pvANFXd1N1KVV0H3ABcEMe4DjhHlpaRm5nD+upNbLFuPmOM2a9IElQpsL/x0R9idaL2KdOXydEjjgTgrY3vpTgaY4xJf5GM4tsIjMMZTt6TsUBcmwUicgUwB2gOWXw18BTwAE6LrR2Yqaoz4nnueFu0/t0uE8e+vHoRF0w4M8VRGWNMeoskQf0F+KWIvKGqLeEr3RF8vwQWxjm2ScC9qnpj2Plm4AzWGAMUAv8Qkc2q+niczx8Xi9a/u1fpjaqmGi6f9xMa25oozh3EReXn2kznxhgTJpIENQM4D3hPRO4HlgI1wEDgaOBHONMc/SrOsR2FMzltuMuAy1W1CqgSkXuAK4G0TFDdlX8HaGhrAmBnQyVzlswFsCRljDEh9pugVLVBRI7DGU7+WyA4nNwD7AKeAH7lJoy4EBEfUA5MFZGZQAPwME6X3zBgRcjmnwJHxOvc8RZJmfeW9hZmv/0oTy1bYK0pY4xxRTSThKrWAVeJyI+BQ3FaTzuB1ara3eSxvTUEp6X2Z5xRhOOBBUDwgeCGkG0bcIbBp6Weyr93x1pTxhizR7RTHbXgtFgSSlW3ASeFLPpQRGYDp7vvc0LW5QK7Ex1TrLor/74vLe0tPLVsgSUoY0y/F+1ksUkhIhNE5LawxX6cch/b6DqjxWF07fJLKyeMmsKVR1/C4NxBeID8zDwyvPuuTLKzoZJv/98Puer5m1i0/t3kBGqMMWkmqhZUElUDPxWRTcAjwETgWuAa4BPgVhFZBuQD19P9YIq0ESz/HrRo/bvM/ehZKhure9wngHX5GWP6t7RsQbkVe8/BGZ1XC8zDGYjxDHALsBwnUS1x1z2UolBjcsKoKTx0zgyOHzl5v9sGu/yMMaa/SdcWFKr6T2Cv3+Cq2oTzwO7VSQ8qzq455nI21W5lfc3mfW4X7PKzZ6aMMf1JWrag+guf18dtX/4pBw0YBoDX0/OPI7TLz+5LGWP6g7RtQfUXuf4cfnHSNdz95u/325KC2Ef5Nbe1sKNhFzvqd7GjvpId9btobGsiEAgQCATA46EoewDFOUUMyi2iOGcgQ/KKycnMjvWfZowxvWIJKg0Mzh3Er069nt8veYK3N76/3+13NlRy5YIbqWqqId+fx+ElQkNLI6sq19LU1kyWz89BhcPIzsiipqmOmuY66ppjG4mf5fMzvGAoA3MKKcweQFH2AAqzC9zve97nZeZapWBjTFxZgkoT2ZnZTDvuv1i2bSVPLVvAmqp9zc3rzOcHsLulnrc3dU1qze0trKnsur/P62Nw7iBK8gYxJLeYwXnF5Ptz8Xo8ePDSHmhn2baVvL91OR2Bji7H+rx6I59Xb9xnPBneDCdxZXVNYM7rQoqyCyjKKWRo3mB8+xlmb4wxYAmqV4KzlO9qqIzbAIby0vGUl45nQ/Vm/rL8eZZsXkaA2CbrKMwq4OaTr2VAVgGF2QX8a8NSnlq2gOXbtdt4n/v05S7JKagoewD/NfliappqqW6qo7qpxmmZNdVS01RHdVMtjW1N7GqoYlfDvme8yvBmMGJAKQcXDufgwhEcXDScMQNHMSDbCjIbY7qyBBWj8FnK4/3M0sFFI7j+Sz/okgQH5RSxax/PToWrba5jVNFBEcfb07yBNU21nbWsetLS1kJ1cx3VjTXUNNe5ycz5CiaxyoYqdjRUsr56E+uru9a/HFZQwmGDv4AMHsNhQ8YwLL/EugyN6ecsQcWou1nKEzFNUfhDvlc9f1PEc/sV5w7qfB1JvD3NGxh6nJ74M/yUZBRTklcMOAnxn2vf2qt12dDayKaarWyo2cyG6i2sr9nM2sr1bK2rYGtdBa99/hYAhdkDOHLoeI4sLaO89DAKswdE9G82xiReIBBgc902lm1bSWVjDReUnU52AgZUWYKKUU+tjUhmL++N7ub283l8eDzQ1tHeuczv83NR+bn7jSt0eXfHDj9OJPbXWhs3+FDGDT60c/u2jnbWVW3k051r0J1r+HTnGmqaanlz/Tu8uf4dAEYPHMmRpWV8sbSMccWHkuGzj64xyVTdWMOy7Z/ysfsVOhPOhJKxTBx2eNzPaf/LY9Sb1kZvBFs74fe+ulsW2vKKJN6ejh1tizDa1mWG18cXig/hC8WHcJacSiAQYFPtVj7atoKPtq1gxY7VfF61kc+rNjJ/5YtkZ2RxeIlwZGkZRw4rozR/SFTxGWP2b3dLPSsqPmNFxSo+rlA21mzpsr4wewBHDD2Mo4YfzhdLJyQkBktQMYpXayMW4d1+oct7Emm8PR07Gr1tXXo8HkYWDmdk4XDOktNoaWthxY7VnQlrU+1Wlm5ZxtItywAYmj+EI0vH88XSMiaUiD27ZUwM6lsaWLljNZ9UrOKTCmV99eYuA7SyfH7KSsZyxNDDKB86npGFwxN+n9gSVIzi1dpIlmTGG8/WZXcjJccP+QLLtq3kw20r+Hj7p2zfvYOXVu/gpdVv4vN4kcFjnNZVaRmHDDwIr8ebkBGXxvRlDa2NfNqZkFbxefVG56F9V4Y3g7HFo5lQMo7DS8Yxtng0mb7MpMboCQ2oPxCRQ4DPX331VQ466KBUh3NACr8HBU5r7cqjL4kqKURynI6ODlZXrnNbVyv5rPLzLv/JCrMKKM0vYU3Vur3u0UUbjzF9VSAQYGdDJbpzLbpzDat2rmVdzaYu/1d8Xh9jBx1CmZuQxhUfij/Dv4+jxsemTZs49dRTAUar6rrQddaCMnGXzHtZXq+3c9DFtw4/i90t9Szfrnzodgfuaqiiprlur2M7x5lvCcockNo62llfvQnduaYzKYWX9/F5vIwpHk1ZyVgOLxHGDT6U7IysFEXcPUtQUbBuosil6l5Wvj+PY0ZO4piRkzoHW/z0H7/qdtudDVU88PZjlJeO5/ASYVBu0T7jidfPP90+R+kWj4lOW0c7m2q28nnVBtZWbeDzqo2sq95IS3trl+3y/LmMKz4UGXwoMngMYwaNSruEFM4SVIQS/WCu2Vtv72UFB1sM7uE4QJeh7MPyS5hQMq6zi6Mop7Bzu3j9/OP5OYpHYrHP9b6lU/Le3VLP3/Q1/rH6dXa31OP3ZVLgz6O2eTetHW17bT+soAQpHtOZkIYPGLrPignpyBJUhJL1YK7ZI14jJXs6zrcmnInX42V5hbJyx2ds3V3B1t0VvLJ2MQAjCkqZUDKOCUPHMfejZ+Py84/X5yheiSXdPtfplBCSnbwDgQA1TbVsr99Jxe5dVNTvZHv9TrbWVbClbvteEz63tLd2ziwzNH8Ihw48mEMHHszogSM5dODB5GflxT3GUMn4WVmCilCqHsztz+J1L2t/xzn7sNNo72hnbdWGziG2n+5Yw+a6bWyu28ZLa97s8djR/vzj9TmKV2JJp891urXmEpG8G1oaqajfSUX9Lrbv3tn5uqJ+Jzvqd+3VLReJ4pyBzD7z9pjiiVWyflaWoCKUqgdz+7t43MuK5Dg+r4+xxaMZWzya88Z/jbb2NlZXrueTCuWTilUsr9Bu98vOyGbp5mWUDRlLrj9nv3HE63MUr8SSTp/rdGvNRXONA4EAu1vqqWysZldDNZWNVSGvq6lscN7Xtzbu85z5/jxK8oopyR9MSd5gSvKKGV5QwvCCUn74/PRup42ubNz3BM2JkKyfVZ9MUCJyJPAQUA6sBa5Q1SWJPGcqH8w1yZfhy+CwIc7EtedPOIPXP3+LPy59aq++/sa2Jn6z+H/xeDwcWnQwE4aOY0LJOA4b/IVuHxiO1+coXoklnT7X+0oIHYEO2trbaOlopbW9jdb2Vlo72mjvaKcjEKAj0NHlqz30dUfHXuudr4C7v/s60N65fXugnZzMHBq6SSh+Xxb3LJ5DXUs9u5t3U9tSz+6WetpDHmPoid+X2Zl4SvIGU5Lvfndf52b2/EdOOv0xkayWd59LUCLiBxYAs4ATgfOBl0RklKrWJuq8fe3BXBNfJ48+Dp83I2Rm+YGceMgUPB4vn1SsYvWuz1lTtZ41Vet57tOX8Xq8jBk0yn3Icc8Q3nh9jiJJLG0d7TS3NdPe0U5boJ2Ojg73+55fpGMGHsyFE87ihVWvUtVUw8DsQs4cdypjBh7Mlrrt+4whEAg4yaKjlZb2Vlrb3e9uEgkua+0Ied3eSktHW7evfV5fl2fVOs8D/Odfro7q+iRSc3sz727+cK/lOZnZFOcMZJBblXpQjlOZ2qlQ7bwvyMqPefaFdPpjIlnJss89qCsiXwUeVdURIcsWA39W1T9GsP8h2IO6Js6a2prRnWuce1jblTVVG7rU1vJ5fXxh0CGdT+XH4yHIV9cs5v+WP091Uy05GdmMGFBKe6CduuZ66lsaaGxr6u0/K61k+jLxezPI9GU6X94MfB4v3uCXN+S1xxuyzoPX69vzusu68G29+Lw+fF4fGV4vm2q2sbxCaWhtJC8zl8nDy5EhY8jz51Dgz6cgK48Cfz75WXn4kzDLQroMIonXw/hw4D2oWwasDFv2KXBECmIxBoDsjKzO6ZUgOI3MGj6pUFZUfMba6g3uQ5Nr+OuKv3eZRuYLg0ZRkj+YvMxcMr0ZZPgy8OChvrWB+hbna2dD1Z4b6rt3sqVue5cHLxvbmlhdua5LTB6Ph2xfFhleH16vjwyPD5/Xi8/jg/39ER/J360e8HszOxOG35dBZuf7jG7X+TvfZ5LhzcDv27Ms05vB8u3Ki6vfoKqphkE5RVww4UxOOuQ/yPBmWH0w4ndPNh5xQOJ7lPpigsoHGsKWNQC5KYjFmG7lZuYwafjhTBrulCDobiLOlTs+Y+WOz2I+R4Y3g2H5Qxg+oJThBUMZXjCUYQUlFGUPIM+fS05mdp977mXc4EP55oTTUx2GiUAykmVfTFD1QPidxFxgdzfbGpMW8vy5TB5RzuQR5QDsbq5nxQ6nlMGm2m1U1O+ksa2ZNveeTUcgQJ4/l/zMXPL8uQzMKexyY314wVCG5Bbj9fatBGRMNPpigloB/CRs2WHA4ymIxZiY5GflMeWgLzLloC+mOhRj0lZfTFCvAR4R+QnwAM4ovnLg2ZRGZYwxJq76XP+AqrYAp+MkpkrgJuA8Vd2R0sCMMcbEVV9sQaGqy4EvpToOY4wxidPnWlDGGGP6B0tQxhhj0pIlKGOMMWnJEpQxxpi01CcHSfSSD2Dbtm2pjsMYY/q9kN/FvvB1/TFBDQO45JJLUh2HMcaYPYYBa0IX9McEtQQ4AdgK7L+AizHGmETy4SSnvWr69blyG8YYY/oHGyRhjDEmLVmCMsYYk5YsQRljjElLlqCMMcakJUtQxhhj0pIlKGOMMWnJEpQxxpi0ZAnKGGNMWuqPM0nETESOBB7CKTG/FrhCVfd6+jmVROQrwF3AWKAC+K2qzhGRLKAOaAnZ/C1V/WoKwuwkIlcAc4DmkMVXA08BDwAX4Mz4MVNVZyQ/wq5E5BKceEPlAK8CZ5Nm11hEpgALVbXEfe9nH9dVRC4E7sR5sv8N4HJVrUhhvCXAfcCpgAf4O3Cdqla56x8HLgTaQg5TrqprUxjzPv+vpeE13h22SQaQBYxQ1S2pvMaWoCLk/sdeAMwCTsQpOf+SiIxS1dqUBucSkZHAPOAynFiPAl4UkXXALqBSVUtTFmD3JgH3quqNoQtFZAYgwBigEPiHiGxW1cdTEGMnVZ0LzA2+F5GJwEvAz4AjSJNrLCIe4HvAPWGrbqOH6yoiZcAjwOnAUuBu4GngyymM92GgBhgNZAJPAA8CF7vrJwHnqeo/Eh1juH3E3OPnIB2vsarmh2yTAbwGvK6qW9zFKbvG1sUXuZOBTFWdpaqtqvo08Anw7dSG1cUhwJOq+qyqdritu9eB43GS1YcpjK0nPcV1GXCHqlap6jqc/1RXJjOw/RGRTJxk9UtV/Yj0usa3AT8Efh22fF/X9VLgeVVdrKpNwHTgeBEZm4p4RcQLdAC3qWq9qlYDfwS+5K7PAQ4jdde8p2u8r89BWl3jbtyA84fArZD6a2wJKnJlwMqwZZ/i/LWUFlR1kar+IPheRAbhTIz7Ac5fQSUiskxEtovI/xOREamK1Y3Ph9NdOlVEtojIahG5UUQG4nR/rAjZPK2utetqoBH4vfs+na7xQ6p6FM5f6QCISBH7vq5loetUtQHYSHKu+17xun9knaeqq0O2Ow/n8wzwRZxupz+KyA4ReV9EzkpCrEF7xeza1+cgra5xKBEZDvwC+IGqdriLU3qNLUFFLh9oCFvWAOSmIJb9EpFC4DngHZzuvnrgXzh9+YLzi/XZlAXoGILzn+XPOF04F+D8hfcjd33o9U6ra+12+f4Mp/UUnHE5ba5xSPdMqGBXTk/XNWWf8R7i7UJErsdJUDe4iwqARTgtg+HAHcBf3HvFCbePmPf1OUjna/wT4B+qGtpaSuk1tntQkavHuRkeKhcIv8GYciIyDicprQAucf8amha2zTRgh4iMVNWNKQgTVd0GnBSy6EMRmY3TPw9dr3e6Xeuv43Q/vRBcoKppd43D1Lvfe7quafkZd7tSZ+MMQvmyqn4KoKov4dz/C5onIt8FzgE+Snqgrn19Dkjfa+zD6f7tUigv1dfYWlCRW4Hz11Cow+jaXZJyInIiTqtpPnCB28+NiNwuIuNDNvW735uSHGInEZkgIreFLfbjxLSNrtc73a71ucBfQrpC0vIah3JHvu3runb5jItILnAwKbzuIlIAvAwcDUwJ/eteRM4WkcvCdgl+flJmP5+DtLvGruPc76+GLkz1NbYWVOReAzwi8hOcYbrn49w/SXU3WScRGQMsBG5S1dlhq8uBySISHP10H/CCqu5IZoxhqoGfisgmnJFNE4FrgWtwBqDcKiLLcLpFrseJOV0cA/xP2LJ0vMbhnqDn6/oksFhETgb+DcwAPlDVVakI1PU0zh/SJ7j3a0L5gPtEZCXwHs6ApeOA7yc3xL30+DkQkXS8xuB8nt8O/YPLldJrbC2oCKlqC07X0/lAJXATztDLdPrlczVOn/EMEdkd8nU3zvDSKmA1sA7nGY2pKYsUUNXNOF0FVwK1OEPkf6WqzwC3AMtxEtUSd91DKQq1O4cA4X36aXeNu9HjdVXVj4Er3Pc7gQnAt1ITJohIOXAGMAWoCPk8b3LjnY/z//ApnM/PT4GzVHVDqmJ29fg5SLdrHOIQ9v48p/waW0VdY4wxaclaUMYYY9KSJShjjDFpyRKUMcaYtGQJyhhjTFqyBGWMMSYtWYIyxhiTluxBXWOiJCKP4UwL05PbcGaRfw0oUNWkTGPjTlfzL+A7sTz4KSIB4GxVXRjBtg8AS1T1z9FHakxkrAVlTPSuw5kVfBhOGRZwHiYNLrsHeMt9Xd/N/olyLfBRL2YlGIYzrVAkbgduF5HiGM9lzH7Zg7rG9IKIHA58DIx26yulKo5sYAPOZKrLk3TOPwEbVPWXyTif6X+si8+YBHDnWuvs4nO7zy7CKVAnOGVGLsUp2TEVZxqZ6ar6hLt/AXAvTgmSAPBPnFLnPZVM+E+gOpicROQQ4HOcqaRmAgcBr+CUM7kHZ2bwrcBV7ozVXbr4ROR1nHLkRwJfxalZ9FtVfTjknH8FHhGRO1S1NeaLZUwPrIvPmOS5C/gxzsScBwPv4ySmo3F+2c8RkWDNpj/gJLKv4ZQkCQAvuiW5u3Mm0F1J7l/hlEf/Ck535DKc7sejcAr/PdzNPkE34HT5TcRJbv8rIqFlzF8Bit1jGRN3lqCMSZ4HVfU1t2TEQpwaQL9QVcVp5eQAo0XkUJwW0cWqusRtFU3FmdDz6z0cezLOBLDh7nCPsRin8NwKVb3fran0IDDSba1153VVfdCNbzpOj0t5cKVbymWte25j4s4SlDHJE1q6vAFYF1KNN1hfJwunLDiABmfwBnYBeexdkyxoKM7s2Ps759qQ96Hn7E7nYAtVrXVfZoZtswso6WF/Y3rF7kEZkzzh92nCa+8EZbjbTsTp2gtV2cM+HYCnF+fsTks3y8LP4QPaozimMRGzFpQx6WclTkslT1VXq+pqnAENvwXG9bDPNmBIkuILNdg9tzFxZwnKmDTj3vN5DnhcRE4QkcOAx3EGV3zaw27v4Yy4SxoRKQRG4RQ+NCbuLEEZk54uwxmKPh8nARQCX1HV6h62fwFntF8yfQmn9fRBks9r+gl7UNeYA4CI5OKUF/+6qr6fpHM+hTMq8FfJOJ/pf6wFZcwBQFUbcO5RXZ2M84nIMJwW24PJOJ/pnyxBGXPg+B1QLiI9DUWPp5uBm1W1p1GFxvSadfEZY4xJS9aCMsYYk5YsQRljjElLlqCMMcakJUtQxhhj0pIlKGOMMWnp/wdlIlzAjVcXvwAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Solution\n",
    "\n",
    "plot(results.I, 'g-', label='simulation')\n",
    "plot(data.insulin, 'go', label='insulin data')\n",
    "\n",
    "decorate(xlabel='Time (min)',\n",
    "         ylabel='Concentration ($\\mu$U/mL)')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Exercise:**  Write an error function that takes a sequence of parameters as an argument, along with the `DataFrame` containing the measurements.  It should make a `System` object with the given parameters, run it, and compute the difference between the results of the simulation and the measured values.  Test your error function by calling it with the parameters from the previous exercise.\n",
    "\n",
    "Hint: As we did in a previous exercise, you might want to drop the errors for times prior to `t=8`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Solution\n",
    "\n",
    "def error_func(params, data):\n",
    "    \"\"\"Computes an array of errors to be minimized.\n",
    "    \n",
    "    params: sequence of parameters\n",
    "    actual: array of values to be matched\n",
    "    \n",
    "    returns: array of errors\n",
    "    \"\"\"\n",
    "    print(params)\n",
    "    \n",
    "    # make a System with the given parameters\n",
    "    system = make_system(params, data)\n",
    "\n",
    "    # solve the ODE\n",
    "    results, details = run_ode_solver(system, slope_func)\n",
    "\n",
    "    # compute the difference between the model\n",
    "    # results and actual data\n",
    "    errors = (results.I - data.insulin).dropna()\n",
    "    return TimeSeries(errors.loc[8:])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "I0       360.000\n",
      "k          0.250\n",
      "gamma      0.004\n",
      "G_T       80.000\n",
      "dtype: float64\n"
     ]
    },
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>values</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>8</th>\n",
       "      <td>11.892725</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>10</th>\n",
       "      <td>-2.299827</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>12</th>\n",
       "      <td>-7.258477</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>14</th>\n",
       "      <td>-7.522652</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>16</th>\n",
       "      <td>-3.263455</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>19</th>\n",
       "      <td>2.001941</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>22</th>\n",
       "      <td>2.652451</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>27</th>\n",
       "      <td>7.047398</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>32</th>\n",
       "      <td>3.742283</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>42</th>\n",
       "      <td>8.929934</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>52</th>\n",
       "      <td>9.543487</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>62</th>\n",
       "      <td>0.764077</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>72</th>\n",
       "      <td>-3.217837</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>82</th>\n",
       "      <td>-10.385020</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>92</th>\n",
       "      <td>-7.384116</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>102</th>\n",
       "      <td>-9.062852</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>122</th>\n",
       "      <td>-3.583691</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>142</th>\n",
       "      <td>-3.586581</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>162</th>\n",
       "      <td>3.171719</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>182</th>\n",
       "      <td>18.704657</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "8      11.892725\n",
       "10     -2.299827\n",
       "12     -7.258477\n",
       "14     -7.522652\n",
       "16     -3.263455\n",
       "19      2.001941\n",
       "22      2.652451\n",
       "27      7.047398\n",
       "32      3.742283\n",
       "42      8.929934\n",
       "52      9.543487\n",
       "62      0.764077\n",
       "72     -3.217837\n",
       "82    -10.385020\n",
       "92     -7.384116\n",
       "102    -9.062852\n",
       "122    -3.583691\n",
       "142    -3.586581\n",
       "162     3.171719\n",
       "182    18.704657\n",
       "dtype: float64"
      ]
     },
     "execution_count": 12,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Solution\n",
    "\n",
    "error_func(params, data)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Exercise:** Use `leastsq` to find the parameters that best fit the data.  Make a `System` object with those parameters, run it, and plot the results along with the measurements."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[3.6e+02 2.5e-01 4.0e-03 8.0e+01]\n",
      "[3.6e+02 2.5e-01 4.0e-03 8.0e+01]\n",
      "[3.6e+02 2.5e-01 4.0e-03 8.0e+01]\n",
      "[3.60000005e+02 2.50000000e-01 4.00000000e-03 8.00000000e+01]\n",
      "[3.60000000e+02 2.50000004e-01 4.00000000e-03 8.00000000e+01]\n",
      "[3.60000000e+02 2.50000000e-01 4.00000006e-03 8.00000000e+01]\n",
      "[3.60000000e+02 2.50000000e-01 4.00000000e-03 8.00000012e+01]\n",
      "[-5.30683629e+02 -6.79500345e-02 -2.29098225e-03  7.86040674e+01]\n",
      "[2.86020929e+02 2.22835658e-01 2.97772381e-03 7.96421423e+01]\n",
      "[2.86020933e+02 2.22835658e-01 2.97772381e-03 7.96421423e+01]\n",
      "[2.86020929e+02 2.22835661e-01 2.97772381e-03 7.96421423e+01]\n",
      "[2.86020929e+02 2.22835658e-01 2.97772385e-03 7.96421423e+01]\n",
      "[2.86020929e+02 2.22835658e-01 2.97772381e-03 7.96421435e+01]\n",
      "[1.40631803e+02 1.53621683e-01 1.89393929e-03 7.90793326e+01]\n",
      "[1.40631805e+02 1.53621683e-01 1.89393929e-03 7.90793326e+01]\n",
      "[1.40631803e+02 1.53621685e-01 1.89393929e-03 7.90793326e+01]\n",
      "[1.40631803e+02 1.53621683e-01 1.89393932e-03 7.90793326e+01]\n",
      "[1.40631803e+02 1.53621683e-01 1.89393929e-03 7.90793338e+01]\n",
      "[1.01576311e+02 8.86883215e-02 9.10961079e-04 7.82514786e+01]\n",
      "[1.01576313e+02 8.86883215e-02 9.10961079e-04 7.82514786e+01]\n",
      "[1.01576311e+02 8.86883229e-02 9.10961079e-04 7.82514786e+01]\n",
      "[1.01576311e+02 8.86883215e-02 9.10961092e-04 7.82514786e+01]\n",
      "[1.01576311e+02 8.86883215e-02 9.10961079e-04 7.82514798e+01]\n",
      "[5.41795211e+01 3.19558519e-02 1.64722970e-06 7.52821022e+01]\n",
      "[9.41838679e+01 7.81822807e-02 7.18446369e-04 7.71896713e+01]\n",
      "[9.41838693e+01 7.81822807e-02 7.18446369e-04 7.71896713e+01]\n",
      "[9.41838679e+01 7.81822818e-02 7.18446369e-04 7.71896713e+01]\n",
      "[9.41838679e+01 7.81822807e-02 7.18446380e-04 7.71896713e+01]\n",
      "[9.41838679e+01 7.81822807e-02 7.18446369e-04 7.71896724e+01]\n",
      "[7.78825680e+01 5.77937352e-02 4.15373341e-04 7.56865719e+01]\n",
      "[7.78825692e+01 5.77937352e-02 4.15373341e-04 7.56865719e+01]\n",
      "[7.78825680e+01 5.77937361e-02 4.15373341e-04 7.56865719e+01]\n",
      "[7.78825680e+01 5.77937352e-02 4.15373347e-04 7.56865719e+01]\n",
      "[7.78825680e+01 5.77937352e-02 4.15373341e-04 7.56865730e+01]\n",
      "[6.35497242e+01 3.52871788e-02 1.00914602e-04 7.17166904e+01]\n",
      "[7.00853351e+01 4.50180592e-02 2.47464050e-04 7.34214050e+01]\n",
      "[7.00853362e+01 4.50180592e-02 2.47464050e-04 7.34214050e+01]\n",
      "[7.00853351e+01 4.50180598e-02 2.47464050e-04 7.34214050e+01]\n",
      "[7.00853351e+01 4.50180592e-02 2.47464054e-04 7.34214050e+01]\n",
      "[7.00853351e+01 4.50180592e-02 2.47464050e-04 7.34214061e+01]\n",
      "[6.37912891e+01 3.44851857e-02 9.62260771e-05 6.83462066e+01]\n",
      "[6.67348714e+01 3.95606416e-02 1.77689353e-04 7.08245914e+01]\n",
      "[6.67348724e+01 3.95606416e-02 1.77689353e-04 7.08245914e+01]\n",
      "[6.67348714e+01 3.95606422e-02 1.77689353e-04 7.08245914e+01]\n",
      "[6.67348714e+01 3.95606416e-02 1.77689356e-04 7.08245914e+01]\n",
      "[6.67348714e+01 3.95606416e-02 1.77689353e-04 7.08245925e+01]\n",
      "[6.34063138e+01 3.36912594e-02 8.92271739e-05 6.53604365e+01]\n",
      "[6.54581258e+01 3.75386018e-02 1.53145128e-04 6.89998360e+01]\n",
      "[6.54581268e+01 3.75386018e-02 1.53145128e-04 6.89998360e+01]\n",
      "[6.54581258e+01 3.75386023e-02 1.53145128e-04 6.89998360e+01]\n",
      "[6.54581258e+01 3.75386018e-02 1.53145131e-04 6.89998360e+01]\n",
      "[6.54581258e+01 3.75386018e-02 1.53145128e-04 6.89998370e+01]\n",
      "[6.37967075e+01 3.44948737e-02 1.05359671e-04 6.52279435e+01]\n",
      "[6.37967084e+01 3.44948737e-02 1.05359671e-04 6.52279435e+01]\n",
      "[6.37967075e+01 3.44948743e-02 1.05359671e-04 6.52279435e+01]\n",
      "[6.37967075e+01 3.44948737e-02 1.05359673e-04 6.52279435e+01]\n",
      "[6.37967075e+01 3.44948737e-02 1.05359671e-04 6.52279444e+01]\n",
      "[6.34409307e+01 3.39762874e-02 1.00577649e-04 6.31260404e+01]\n",
      "[6.34409316e+01 3.39762874e-02 1.00577649e-04 6.31260404e+01]\n",
      "[6.34409307e+01 3.39762879e-02 1.00577649e-04 6.31260404e+01]\n",
      "[6.34409307e+01 3.39762874e-02 1.00577650e-04 6.31260404e+01]\n",
      "[6.34409307e+01 3.39762874e-02 1.00577649e-04 6.31260413e+01]\n",
      "[6.28791031e+01 3.28824014e-02 8.29220012e-05 6.00474015e+01]\n",
      "[6.30778597e+01 3.32816079e-02 8.95761488e-05 6.11905313e+01]\n",
      "[6.30778606e+01 3.32816079e-02 8.95761488e-05 6.11905313e+01]\n",
      "[6.30778597e+01 3.32816084e-02 8.95761488e-05 6.11905313e+01]\n",
      "[6.30778597e+01 3.32816079e-02 8.95761502e-05 6.11905313e+01]\n",
      "[6.30778597e+01 3.32816079e-02 8.95761488e-05 6.11905322e+01]\n",
      "[6.28214056e+01 3.28104757e-02 8.25234881e-05 5.95114505e+01]\n",
      "[6.28214065e+01 3.28104757e-02 8.25234881e-05 5.95114505e+01]\n",
      "[6.28214056e+01 3.28104762e-02 8.25234881e-05 5.95114505e+01]\n",
      "[6.28214056e+01 3.28104757e-02 8.25234894e-05 5.95114505e+01]\n",
      "[6.28214056e+01 3.28104757e-02 8.25234881e-05 5.95114514e+01]\n",
      "[6.27778545e+01 3.27546603e-02 8.21350019e-05 5.93234986e+01]\n",
      "[6.27778555e+01 3.27546603e-02 8.21350019e-05 5.93234986e+01]\n",
      "[6.27778545e+01 3.27546608e-02 8.21350019e-05 5.93234986e+01]\n",
      "[6.27778545e+01 3.27546603e-02 8.21350032e-05 5.93234986e+01]\n",
      "[6.27778545e+01 3.27546603e-02 8.21350019e-05 5.93234995e+01]\n",
      "[6.27660608e+01 3.27361764e-02 8.19166237e-05 5.92757881e+01]\n",
      "[6.27660617e+01 3.27361764e-02 8.19166237e-05 5.92757881e+01]\n",
      "[6.27660608e+01 3.27361769e-02 8.19166237e-05 5.92757881e+01]\n",
      "[6.27660608e+01 3.27361764e-02 8.19166250e-05 5.92757881e+01]\n",
      "[6.27660608e+01 3.27361764e-02 8.19166237e-05 5.92757890e+01]\n",
      "[6.27629251e+01 3.27314482e-02 8.18644601e-05 5.92648489e+01]\n",
      "[6.27629260e+01 3.27314482e-02 8.18644601e-05 5.92648489e+01]\n",
      "[6.27629251e+01 3.27314487e-02 8.18644601e-05 5.92648489e+01]\n",
      "[6.27629251e+01 3.27314482e-02 8.18644614e-05 5.92648489e+01]\n",
      "[6.27629251e+01 3.27314482e-02 8.18644601e-05 5.92648498e+01]\n",
      "[6.27620378e+01 3.27300853e-02 8.18488907e-05 5.92615370e+01]\n",
      "Both actual and predicted relative reductions in the sum of squares\n",
      "  are at most 0.000000\n"
     ]
    }
   ],
   "source": [
    "# Solution\n",
    "\n",
    "best_params, details = leastsq(error_func, params, data)\n",
    "print(details.mesg)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>values</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>I0</th>\n",
       "      <td>62.762</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>k</th>\n",
       "      <td>0.0327301</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>gamma</th>\n",
       "      <td>8.18489e-05</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>G_T</th>\n",
       "      <td>59.2615</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>G</th>\n",
       "      <td>&lt;function interpolate.&lt;locals&gt;.wrapper at 0x7f...</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>init</th>\n",
       "      <td>I    62.762038\n",
       "dtype: float64</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>t_0</th>\n",
       "      <td>0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>t_end</th>\n",
       "      <td>182</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>dt</th>\n",
       "      <td>1</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "I0                                                  62.762\n",
       "k                                                0.0327301\n",
       "gamma                                          8.18489e-05\n",
       "G_T                                                59.2615\n",
       "G        <function interpolate.<locals>.wrapper at 0x7f...\n",
       "init                         I    62.762038\n",
       "dtype: float64\n",
       "t_0                                                      0\n",
       "t_end                                                  182\n",
       "dt                                                       1\n",
       "dtype: object"
      ]
     },
     "execution_count": 14,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Solution\n",
    "\n",
    "system = make_system(best_params, data)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>values</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>success</th>\n",
       "      <td>True</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>message</th>\n",
       "      <td>The solver successfully reached the end of the...</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "success                                                 True\n",
       "message    The solver successfully reached the end of the...\n",
       "dtype: object"
      ]
     },
     "execution_count": 15,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Solution\n",
    "\n",
    "results, details = run_ode_solver(system, slope_func, t_eval=data.index)\n",
    "details"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAagAAAEYCAYAAAAJeGK1AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjAsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+17YcXAAAgAElEQVR4nO3deXxU1f3/8ddkg4QAgbAvAlL8AEoUkBQXwAW/rohba5G6fv2WqlVba2ut1rW4fEWLIlYsLl8Vsf2JimJdERXrAgIKCBxEZA0QIAkEsifz++NO4iQkIRMmmUnyfj4eeUzm3DtzP9wM+eSc+7nn+Px+PyIiItEmJtIBiIiIVEcJSkREopISlIiIRCUlKBERiUpxkQ6gsZlZK2AEsA0ojXA4IiItXSzQHVjsnCsM3tDiEhRecloY6SBERKSSUcCnwQ0tMUFtA5g1axbdunWLdCwiIi3a9u3bmThxIgR+NwdriQmqFKBbt2706tUr0rGIiIjngEsuKpIQEZGopAQlIiJRSQlKRESikhKUiIhEpZZYJCEiTcjevXvJzMykuLg40qFIPcTHx9OlSxfatWsX8muVoMJk4cZFzF4+l915WaQmdWRC2nhG9UmPdFgiTdrevXvZsWMHPXv2JDExEZ/PF+mQJAR+v5/8/Hy2bt0KEHKS0hBfGCzcuIgZi2exKy8LP7ArL4sZi2excOOiSIcm0qRlZmbSs2dPkpKSlJyaIJ/PR1JSEj179iQzMzPk1ytBhcHs5XMpKi2q1FZUWsTs5XMjFJFI81BcXExiYmKkw5BDlJiYWK8hWiWoMNidlxVSu4jUnXpOTV99f4ZKUGGQmtQxpHYRkca0efPmSIdQL0pQYTAhbTwJsQmV2hJiE5iQNj5CEYlItHjyySe56aabwvqeW7ZswczYu3fvQfedNWsWDzzwQMXzoUOH4pwLazwNRVV8YVBeracqPhGp6te//nVEj5+VlYXf7694vmzZsghGExolqDAZ1SddCUmkBSsrK+OBBx5g3rx5+P1+Bg8ezJ133sncuXNZvXo1TzzxBNOmTSMjI4M9e/bw+eef06tXL+69916eeeYZFi5cSK9evZgyZQpmxrRp0ypeB16v6dRTT2Xx4sUHHPvdd9/lqaeeYvPmzfj9fsaOHcs999zDhx9+yIwZMygrK+Pcc8/ljTfewMx4/fXXGTRoECtXruTBBx9k9erVpKamctlll5XPLM6ll17KsGHD+OSTT9i4cSMDBw7k3nvvpX///o12TpWgRKRJuf+T6SzbtrJRjjW0+1HcOvq6Ou37/vvv88knn/D222+TlJTEHXfcwYwZMw5Y1ueNN95g+vTpTJs2jeuuu45f/vKXTJ8+nSlTpvCnP/2Jxx9/nGnTptU5xq1bt/LHP/6Rp59+mmOPPZaNGzdy8cUX88EHH3DmmWeydu3aSomuXFZWFldccQXXXXcdzzzzDGvXrmXSpEm0b9+ec845B4DXX3+d5557js6dO3PjjTcybdo0pk6dWufYDpWuQYmIhEHbtm3JzMxk7ty5bN++ncmTJzN58uQD9ktLS+Okk04iNjaW9PR0+vfvz5gxY0hISOD4449ny5YtIR23c+fOzJs3j2OPPZbc3FyysrLo0KHDQe87mj9/Pp07d+bKK68kPj6eI488kssuu4w5c+ZU7HPuuefSr18/kpOTOf3009m0aVNIsR0q9aBEpEmpa4+msR1//PHccccdvPzyyzzwwAP07t2bW2655YD9UlJSKr6PiYmpNLtCTEwMZWVlIR03Pj6eOXPm8Morr9C6dWsGDx5MYWFhpetO1cnKyqJHjx6V2nr27Mm2bT+uG5iamlrxfVxcHKWlByzZ1KCUoEREwmDz5s0MHjyYl19+mdzcXF566SV++9vfcuWVV1bar673BMXExFS6uTUnJ6fa/d566y3efPNN5syZQ9euXQGv53Mw3bt3JyMj44B/Q6dOneoUX2PQEJ+ISBh88cUX/OY3vyEjI4Pk5GTat29P27ZtiY2Nrdf79evXj2XLlrF+/Xry8vJ49tlnq90vNzeX2NhYEhISKC4u5oUXXsA5V5HcEhISyM3NPeB1Y8aMITs7m+eee47i4mJWrVrFCy+8wLhx4+oVb0NQghIRCYMLL7yQsWPH8vOf/5xhw4bxr3/9i8cee6zesyiMHTuWs88+mwkTJnDWWWdx3HHHVbvf+eefz+DBgxk7diyjR4/miy++4JxzzuG7774D4KSTTmLjxo2MGTOm0uvat2/PzJkz+eCDDxg5ciS/+c1vuPrqq7n44ovrFW9D8B1snLK5MbO+wA/z58+nV69ekQ5HRGqxevVqBg0aFOkwJAxq+lmWl88D/ZxzG4K3qQclIiJRKSqKJMwsHZjnnOsSeN4FeBQ4FfABbwM3OueyA9t/DtwHdAc+Bq5wzoU+l7uIiEStiPagzMxnZlcD7wHBk9nNBEqAfsAAoAMwPfCawcDTwBVAKvAd8HLjRS0iIo0h0kN8dwPXAH8tbzCzGKAMuNs5t985lwP8AzgxsMsvgTedc5865wqAW4ETzGxA44YuIiINKdIJ6knn3HDgq/IG51yZc+4859y6oP3OA8pnOBwMrAraPw/YDAxphHhFRKSRRPQalHMu42D7mNnNeAnq+EBTMpBXZbc8ICm80YmISCRFRZFEdcwsHpgGjANOcc6tCWzaD1RdAzoJ2NeI4YmISAOLygRlZm2BN4G2QLpzbmvQ5lWABe2bBBxG0LCfiIg0fVGZoPCq8mKAUYFrTMFeAj41s5OAz4H7gWXOubWNG6KISMtRWlpKZmYm3bt3b7RjRrpI4gBmlgacBaQDmWa2L/C1BcA5twK4CngS2AUcCfwsUvGKiGRkZDB06NBq57wLl0svvZTnnnsOgKuvvppZs2Yd8nu++uqrjB8/vk773nTTTbz77ruHfMxQREUPyjn3EZAS+H453s25te0/B5hT2z4iIgALNy5i9vK57M7LIjWpIxPSxod99esePXo06lLqM2fObLRjlcvKymr0Y0ZdD0pEJFwWblzEjMWz2JWXhR/YlZfFjMWzWLhxUViPs2XLFsyMvXv3smXLFoYOHcqzzz7LiSeeyHHHHcedd95Zsc7TRx99xNlnn82xxx7LuHHjeP311w94j3LBvaZgwe2XXnopf/vb3zj//PMZNmwYl1xyCd9//321cebk5HD99dczbNgwzjjjDFau/HFlYr/fz+OPP86ZZ57J0KFDGT16dMUxJk+ezFdffcWUKVO45557AHjppZcYN24cw4cP57jjjuOhhx461NN4gKjoQYmINITZy+dSVFpUqa2otIjZy+eGvRcVLC8vD+ccH3zwAevXr2fixImccsopjBo1iptvvpmHH36YMWPG8Nlnn3HDDTdw2mmnHdLx6ro0+x133EFRURGffPIJu3fv5qqrriI5ORmAefPmMXfuXF588UW6dOnChx9+yPXXX89ZZ53Fbbfdxpo1azj11FO54oorWLp0KVOnTmX27Nn079+f5cuXc8kll3D66aeTlpZ2SP+WYOpBiUiztTuv+mGpmtrD6Ve/+lXFCrdmxqZNm4iJiaFNmza89dZbfPXVV6Snp7No0SLatGlzSMeqy9LshYWFFUknOTmZPn36cOmll1ZsP/nkk5k1axZdu3Zl165dxMfHU1paWu3Q3qBBg3j99dfp378/2dnZFBQU0KZNm4MuMx8q9aBEpNlKTerIrmqSUWpSxwY/dvDKtOW/7AGeeeYZpk2bxrXXXktpaSkXXnghf/jDHw7pWHVZmj0nJ4fi4mK6detW0Ra85FBJSQn3338/n332GV26dKnoCVW3JFNsbCwzZszg3XffpUOHDgwePDjkperrQglKRJqtCWnjmbF4VqVhvoTYBCak1a1yLdzy8/PJzMxk6tSplJWVsXTpUm644QaOPPJI0tO9Ice6LPNeHx06dCAhIYGMjIyK5Lljx46K7Y888giFhYV8/PHHtG7dmj179vDKK69U+17PPvssq1at4r333qNdu3b4/X5GjBgRtljLaYhPRJqtUX3SmTRiIp2SOuIDOiV1ZNKIiQ16/ak2paWlXHPNNbz11lv4fD66deuGz+cjJSWF1NRU2rZty+uvv05paSlvv/12jcUO9ZGQkMA555zD1KlT2bNnD1u2bOH555+v2J6bm0urVq2IjY1lz5493HfffQCVlo7ft29fxb7x8fHExcWRn5/PI488Qm5uLkVFRQce+BCoByUizdqoPukRS0hVJScn89hjjzFlyhRuv/12kpOTmThxYsVy7Pfddx8PP/ww06dP5+STTz7k4omqbr/9du666y5OPvlkUlJSGDt2LF9++SUAN954I3/6059IT0+nbdu2nHXWWZgZa9euJS0tjXHjxnHPPffwww8/cNttt7F69WpOOOEEkpKSGD16NCeccELFMvPhoiXfRSRqacn35kNLvouISLMR8hBfYEHBgUAXoBTYDqxzzrWsrpiIiDSoOicoMxsN3AiMxZtlvJwfyDazd4AnnHOfhTdEERFpiQ6aoAJLqc/AW9LiNeACvKUtduMNEXYGjgZGAy+b2ffAJM0uLiIih6IuPagXgXucc2/VsH1z4Guemd2Ct/rti3izkYuIHBK/34/PV+v80RLl6luMV5cENbKu15cC+71mZq/XKxoRkSDx8fHk5+eTlJQU6VDkEOTn5xMfHx/y6w5axVeX5GRmPczsslBeIyJyMF26dGHr1q3k5eXV+69wiRy/309eXh5bt26lS5cuIb8+XDfqDgGeBZ4/2I4iInXVrl07wFsQMHgKIGk64uPj6dq1a8XPMhSaSUJEolq7du3q9ctNmj7dqCsiIlFJCUpERKJSXe6DOqsO7zM8DLGIiIhUqMs1qHl1fC+V2IiISNgcNEE55zQMKCIija4+k8UmAf2AVlU2+Z1zy8ISlYiItHghJSgz+yXwJJAIVJ17xA/E1icIM0sH5jnnugSeJwCPAxfhzZj+iHPu/qD9fw7cB3QHPgaucM5l1ufYIiISnUIdvrsfeBo4HC85BH/1CPXgZuYzs6uB94CEoE13Awb0B0YAl5fPVGFmgwMxXAGkAt8BL4d6bBERiW6hDvG1Ax53zm0M0/HvBs4G/grcHtR+OV6vKBtvKY8pwCS8mSp+CbzpnPsUwMxuDewzwDkX3vWGRUQkYkLtQb2A13MJlyedc8OBr8obzCwFr0e2Kmi/NXjTKQEMDt7mnMvDm019CCIi0myE2oN6CFhqZhOBDUBZ8Ebn3CmhvJlzLqOa5uTAY15QWx6QFLQ9j8qCt4uISDMQaoJ6AdgHvMWBSSJc9gceE4PakgLHLd+eSGXB20VEpBkINUGNAH7qnFveEMEAOOeyzWw7XpHE1kDzQH4c1lsV2AZUlL0fRuUhQRERaeJCTVAOSGmIQKp4AbjTzJbjDendDDwa2PYS8KmZnQR8jldZuExLzIuINC+hJqj7gefM7HHge6DSAi3OuX+HKa47gIeBb/EKOZ7Cu/8K59wKM7sq8Lwn8CXwszAdV0REokSoCWp24HFKNdvqfaOuc+4jgnpmzrkC4LrAV3X7zwHm1OdYIiLSNISUoDQvn4iINJa6LLfxHd5MD+8D851zuQ0elYiItHh16UGdDZwGXAk8Y2ar+TFhfeGcK6vtxSIiIvVRl+U21gJrgelmFgsch5ewHgYGmNlCAgnLObeuIYMVEZGWI9RrUKXAp4GvO82sPXAqXsK6CRgQ9ghFRKRFqss1qME1bCoGcpxzrwKvhjUqERFp8erSg1qJV0Jedf0nAL+ZbQMecM49HtbIRESkRatLgupXQ3sM0AE4AbjbzPY5554LV2AiItKy1aVIora1n37Am908H7gReC5McYmISAsXrhtvP0EFEiIiEkbhSlAxQGGY3ktERCRsCeoKYEmY3ktERKROZebX1rApBmgPHI93L9SpYYxLRERauLpU8f2hhvZiIBtYBoxwzq0IW1QiItLi1aWKr6YycxERkQYT0lRHZvZrYCmwwjmXH9SeCBDcJiIicihCXbDwVqA3UGpmDm94bylQAtyASs0rLNy4iNnL57I7L4vUpI5MSBvPqD7pkQ5LRKTJCHWy2D5m1hEYDgwDRgB3AcnAhnAH11Qt3LiIGYtnUVRaBMCuvCxmLJ4FoCQlIlJHofagcM5l4a0F9T6AmaXgTRb79/CG1nTNXj63IjmVKyotYvbyuUpQIiJ1dMj3QTnncoC/APcdejjNw+68rJDaRUTkQCElKDPrVcOmnUBN21qc1KSOIbWLiMiBQh3i22Rmu/BmjVga+NoKXAe8EubYmqwJaeMrXYMCSIhNYELa+AhGJSLStISaoNKAoYGv44Fr8WaTAPiPmT0GfA187ZxbGrYom5jy60yq4hMRqb9Qq/hW4i1g+EJ5m5n148ekNRS4EOgGxIYvzKZnVJ90JSQRkUNQl7n4UgKFENVyzv2Aty7Uq0Gv6XKogZnZSOAxwPCucT3gnJtpZgnA48BFQCnwiHPu/kM9noiIRJe6FEl8bGa3mFn7g+1oZp3M7DZg/qEEZWYxwFzgMedce2AC8LiZHQ3cjZe0+uPdh3W5mV12KMcTEZHoU5chvhOAvwJbzOw/wDvAt8AuwAd0Bo4GxgCjgP8LvOZQdAC6AD4z8wF+vNkqioDLgSucc9lAtplNASYBzx/iMUVEJIrUZbLYfcBvzewBvERwCd61pvJrTMV4Ux69BfyPcy7jUINyzu02s8fxkt2zgWPdCGwDugOrgnZfAww51GOKiEh0qXORhHNuO97w2t2BIbhUoMw5tzvcQQXevwAvGc7Bqxh8FSi/FpYXtHsekBTuGEREJLJCnuoIwDlXhle40FAuAE5wzpWvRfWxmT2NN7wHkBi0bxKwrwFjERGRCAjXku/h1htoVaWtBC8pbscrkig3kMpDfiIi0gzUqwfVCN4D7jezXwH/wJs5/X+Aq4FNwJ1mthxvFvWbgUcjFaiIiDSMqOxBOee+xRvmm4R33ekl4E/OubnAHXg3C38LLMa7RvVkhEIVEZEGUq8elJnFB17rC253zuVV/4rQOef+Dfy7mvYCvLn/rgvXsUREJPqEuuT7SGAGcFQNu7To6Y1ERCR8Qu1BTQX2AOcBe8MfjoiIiCfUBDUEGOmcW9EQwYiIiJQLtUhiNdCjIQIREREJFmoPahrwDzObBnyHNzdehUBhg4iIyCELNUE9G3h8sJptflQkISIiYRLqgoVRed+UiIg0P/W9D+pU4Ei8a1irgfnOuZJwBiYiIi1bqPdBdQNex5t6aAPejbp9gDVmNtY5lxn2CEVEpEUKdcjuUbxJW/s5545wzg0A+gJZwCNhjk1ERFqwUBPUGcANzrmt5Q2BBQp/D5wVzsBERKRlCzVBFeBV61WlCj4REQmrUBPUe8AjZta1vCHw/cPAu+EMTEREWrZQq/j+AHwIbDSzjYG2PsByvOXZRUREwiLU+6C2m1kacDpemXk+sNo590FDBCciIi3XQROUmZ0FvO+cKw58D941p5WB7xPK2zXVkYiIhEtdelDzgG5AZuD7mqhQQkREwuagCSp4eiNNdSQiIo0lpIRjZh+aWUo17Z3NbEn4whIRkZauLtegTgIGB56OASaZWW6V3QYB/cMbmoiItGR1uQa1G7gZb949H3AdUBq03Q/sw5tNQkREJCzqcg1qBXA4gJktAC5wzmU3dGDR6tVVb7O3cB+XH3MRPp8v0uGIiDRbod4HdXJ17WaWAAx3zn0elqii2MKNi9i6dzvW6XCO6z080uGIiDRboS63MRKYgXdNqmqBhT/U9zvIsboDfwdOxpsD8Cnn3F8CyfBx4CK8ocZHnHP3h+u4B3P2Eafw1FcvMXv5XEb0PIa4GFXWi4g0hPost7ET+AXeLBK/BG7FuwY1IbyhMRfYBnQFRgKXm9klwN2A4RVljAi0XxbmY9fo5H7H071tF7bv28n87z9trMOKiLQ4oSaoNOAm59wcYBmwwzn3v8CNwG/DFZSZ/RTvutcNzrkC59wPwEnAAuByYLJzLts5twGYAkwK17EPJjYmlkvSzgPglW/foqC4oLEOLSLSooSaoEqAvYHv1wLHBL5fgDc3X7gMB1YAd5nZVjP7Hjgfr9fWHVgVtO8aYEgYj31Q6T2PYUBqP/YU5vKmpiEUEWkQoSaoL4FrzSwG+AY4M9B+FFAUxrg6AqOAYrye1AV4pe7nBrbnBe2bBySF8dgH5fP5mJh2PgBvug/IKdh7kFeIiEioQk1QtwJXAjcBLwCDzWw98C/gpTDGVQjsdc7d5ZwrdM59A8zEG94DSAzaNwnvGlijGtxlAMN6DKGgpJD/t7K2KQpFRKQ+QkpQzrnFQF/g/wL3Qg3HuwZ0FfC7MMa1BkgKVOyViwOyge14RRLlBlJ5yK/RTEw7jxhfDB+s/5QN2ZsjEYKISLMV6lx884DDnHM7wVsfyjn3hHPun8656paCr6/38aoFHzazBDMbAvw3MBuv53anmXUys754Q38vhPHYdda7fQ/OGHASfr+fp5f+E78/nKdARKRlC3WIbyTedaEG5ZwrwJv373C8UvN3gP8NVA/egbcW1bfAYmAO8GRDx1STnx95Du1btcXt+p6FGxdFKgwRkWYn1Btr/wY8b2Z/A9bjVdVVcM6FbajNObceOLua9gK8+QCvC9exDkVSQiITjz6fJxY9z6xvXmNEz6NJjG8d6bBERJq8UHtQ9+L1ov6J13tZiVcOXv7YIo3u+1MGpPYju2APc1ZpUWERkXAINUH1q+br8KDHFinGF8NVwy7Gh4+33Hw278mIdEgiIk1eqAnqWWCPc25j8BfevUivhj+8pqN/xz6c2v9ESv1lzFg8izJ/WaRDEhFp0uqzYOGvzKzqfUdasBD4Zdr5LNm6nLW71/Peuk84Y8BJkQ5JRKTJqs+Chb9BCxZWKykhkf8e/gum/GcGLy1/nWN7ptEpqWOkwxIRaZK0YGGYpfc6hvSex7Bo69fMXPIyt5x4jRY2FBGph3otWGhm8YHX+qpsz6vudS3NVcMvZkXmGpZmrODzzUso9Zcxe/lcdudlkZrUkQlp4xnVJz3SYYqIRLVQFyz8KfAU3uSwwXx4Q31avQ/omJjCpUdfwFNfvcTfF72InzKKSr37m3flZTFj8SwAJSkRkVqEeqPuo8Ae4Dx+XHZDqnHK4Sfw+ealrNix5oBtRaVFzF4+VwlKRKQWoSaoIcDIwHUpqUWML4Zr0y/jmjf/XO32XXlZXPzPazTkJyJSg1Dvg1oN9GiIQJqj1KQOtE1oU+N2Pz8O+WkePxGRykLtQU0D/mFm04DvqLJIoXNO8/xUcfnQn/HEov+jrJaZzjXkJyJyoFAT1LOBxwer2aYiiWqM7vtTCksKeXrJy5RRc5LanZfViFGJiES/UMvMQx0SFOC0n4ymW9su/PWjx/DXkKRSdUOviEglofagADCzwXir2r4HdAE2hHnBwmZnSNeBXDD4zGpnO0+ITWBC2vgIRCUiEr1CXVG3nZn9G295jVeArnhrRH1jZr0aIL5m5WdHns3gzgMAiI+Jwwd0SurIpBETdf1JRKSKUIfsHgZaAb3wZjAHuAHIAaaGMa5mKSYmhhuOu4p2rZIpLivhZ0eNY0LaeGYvn8vF/7yGa9+8TdV8IiIBoSaos4E/OOcqFjxyzm0CrgdOCWdgzVXHxBSuH3klPnz8a+Wb/H3RC+zKy1LJuYhIFaEmqGSqLPMeEFuP92qxju42mJ8ddQ4AJWUllbaVl5yLiLR0oSaVd4C7ApPFAvjNrDMwBXg/rJE1cxcMPqPGbSo5FxEJPUFdD/TGWyMqCfgA2Ay0B24Mb2jNW4wvhtTEDtVuU8m5iEiICco5t8M5dzxwLl5xxGPAOODY4OtSUjeXHH0e8THxldpUci4i4gl1uY044A5gq3NueqDtS2CkmU12zpU1QIzNVnlp+fPLXmFPYS4AJ/QerpJzERFCv1F3CjAe+FVQ23TgTqA1cFuY4moxRvVJZ1SfdD5c/x+eXPwiCzZ8zpJtK8gt3KeZzkWkRQv1GtTFwCXOuYqCCOfc88DlwBVhjKuCmaWY2SYzuyLwPMHMnjKzLDPbaWa3NsRxG9sph5/AsT3SANhbuE9l5yLS4oWaoJKAfdW0Z+MVSjSEJ4GeQc/vxptmqT8wArjczC5roGM3qg05Ww5oU9m5iLRUoSaoBcD/mllFmZmZpQCTgY/CGFf5e18OtAOCF0i8HJjsnMt2zm3AG3acFO5jR0JN5eUqOxeRlijUa1A34JWWbzWzjXhLbPQB1uFdmwobM+uHd23reLz7r8qTYXdgVdCua/BW+m3yUpM6squaZNQhMSUC0YiIRFaoZeabgKOAi4CZwN+BC4BjnHM/hCsoM4sFXgRuds5tD9qUHHjMC2rLwxt6bPImpI0nITah2m37Cvc3cjQiIpEV8vREzrkivGU2puMlqU+A1mYWziTxF+9Q7tUq7eW/pROD2mq6LtbkjOqTzqQRE+mU1BEf3rx9Ka3akZWfw70fPUpuYbP4Z4qI1Emo90GNBGbg9aKC+Qjvirq/AHqY2QWB522BJ4B0YDtekcTWwLaBVB7ya9LKy87LZeXlcNeCR/ghZzN3L5jKX066gfat20UwQhGRxhHqNaipwB7gPGBv+MPxOOcGBj83s6+Bqc6558xsH3CnmS3HG/K7GXi0oWKJtI5JKdx1yk3cs2Aqm/Zs5e4FU7njpBtJSWyookkRkegQaoIaAox0zq046J4N5w68dam+xRuifAqvFL3Z6pjoJal7F0xl895t3LngEW4fcwOd26RGOjQRkQYTaoJaDfSgctl3g3POHRP0fQFwXeCr2Vu4cRGzl89ld14WHRJTSE3swLbcTG6f/xC3jb6ew1J6HvxNRESaoFAT1DTgH2Y2DfgOKAre6Jz7d7gCEy85zVg8i6JS7zRn5eeQEBtPj7ZdycjdwR0fPswfT7yGwV0GRDhSEZHwCzVBPRt4fLCabeEskhBg9vK5FcmpXFFpMYWlRaT3OoZFW75m8seP8ZuRV3Bc7+ERilJEpGGElKCcc1o1txHVNINEVl4208/+K88s+yfvrfuEv302k61HbefCwWfh8/kaOUoRkYYRag8KADM7FTgSr0hhNfAho94AABlASURBVDDfOVdS+6skVDXNLJGa1JGYmBj+e9gv6NqmMy9+8yr/WjmPzXu2cW36ZbSKq/5mXxGRpiSkHpGZdTOzL4C3gd/gFSq8CSw1sy4NEF+LVt3MEsELGvp8PsYNHMsto64lMa41n29ewp0fPkxWXk4kwhURCatQh+weBUqAfs65I5xzA4C+QBbwSJhja/GqzizRKakjk0ZMPGB9qGE9jmLy2D/StU0n1mdv4tb3H2Dd7g0RiVlEJFxCHeI7AzjZOVc+iwPOuQwz+z3wfs0vk/qqOrNETXq1787k027hkf88xaqd33HHhw9zxdCfcVr/UbouJSJNUqgJqgCvWq8qVfA1ouB7o8pX3QWv6m9XXhat41pRUFLIzCWzWbXzOyYdO5HE+Nb1fm+t6CsikRBqgnoPeMTMfuGc2wFgZl3xZnZ4N9zByYGq3hu1Ky+LJ758Hp8PSspKASgoKSQuJhb88Nmmr/ghexO/O+5q+nboHfJ7z1g8C0BJSkQaXajXoP4AdAU2mpkzMwdsBNoAN4Y7ODlQdfdGlfpLK5JTuZKyUtq2Suaw9j3ZlpvJrR88yNzV71FWVhbSe2tFXxGJlFDXg9oOpOGtB/U08DhwjnMu3Tm3rQHikypCWV03p2APk8f+kf/6yWhKy0qZtfw17v5oKjv37w7pvbWir4hEQp0SlJm1MrNfm1kH51yJc26ec+5/gVbeZtONN40kNaljSPu2ikvg6uETuHX0daS0bsfqnd9x8zt/5b11H1PmLztg/0M9pohIuBw0QZlZe7xFCf+Gtw5TsO540x59aGZtwx+eVFXdvVGxvljvmlOQ4PulAIZ2P4qLjjybhNh48ksKmLnkZW56+24y9v64YPHB7rsSEWlMdelB/QVv1doBzrkvgjc4536Pt3hhZ+DP4Q9Pqqru3qhrf3oZ16RfVuv9Ugs3LuL5r+dQVFpc0ZaRm8nv3r6HV1e9TUlpSZ3vuxIRaQx1qeK7AJjknNtS3Ubn3AYzuwV4CLg1nMFFu0iVZNd0b1Rtx66uAALAj5+XV7zBZ5uW8D/HTqjzfVciIg2tLj2obsDag+zzNd46US1GeUn2rrws/PxYkr1w46JIh1at2goduiZ3ZtOerfxl/hSmffEs2fl7GjEyEZHq1SVBbQaOOMg+A4DtB9mnWWlqJdk1FTp0SurIlNNv58LBZxEfE8fCjYu48d938saa9ygp1fy/IhI5dUlQ/wLuqqlSL9B+FzAvjHFFvaZWkl1bAUSruAQuHjKOR868g2N7Hk1BSSEvfvMav3/3XpZkrMDvr27yEBGRhlWXa1D3A+cBS8zsMeArYA/QARgBXI83zdG9DRVkNKptKYxoVH5dqbZrZl2TO/PHE3/N19tW8dyyf5GRu4MHFz7BoM4/4ZK087BO/SMVvoi0QL66/HUcKCF/ELgEKC8n9wG7gReAe51z2Q0VZDiZWV/gh/nz59OrV696v0/VaYHA65E0l6q3ktIS3ln3Ma+tepvcov0AHNsjjV8MOZfDUnpGODoRaS62bNnCqaeeCt4qGRuCt9VpLj7nXC5wrZn9Fjgcr/e0C1jnnGuR4z916ZE0ZXGxcZxjp3JKv+N5033AvLXz+SpjOUsyVjCqTzoXDD6DHu261fh6TTorIoeqTj2o5iRcPaiWJqdgL6+uepv3v19IaVkpPnz8tPdQLhh0xgGT0Db33qWIhE9tPahQJ4uVFiqldTuuGnYxj555F2MPP5HYmFi+2LyUP753Hw98Mp21u9ZX7NvUKhxFJDqFutxGozGz04AH8ErYM4GHnHMzAlWDj+NNWFsKPOKcuz9ykbYsXZI78asRE7noyLN5w73PB98vZOm2lSzdthLr1J+zjji52uIRiN4KRxGJTlGZoMysNzAHuByYCwwH3jWzDcBJeHMC9gfaA++Y2Vbn3PMRCbaF6piUwhVDf8YFg87grbUf8u66j3G7vsft+p4YX8wBE9FC9FY4ikh0isoEBfQFXnLOvRZ4vtjMPgJOwEtaVwSqBrPNbAowCVCCioB2rdsyIW085w86nY82fMHb3y1gW27mAftp0lkRCVVUJijn3EJgYflzM+sIjMIrae8OrArafQ0wpFEDlAO0jm/NGQNO4r9+Mpqvt61i1jevsXlvRsX2dq2SKSguJK84n6T4xAhGKiJNRVQmqGCB5T7eAL4ElgSa84J2ycObbV2iQIwvhmE9jmJYj6PYsncb87//Dx9v+IJdeVn8Y8lLPP/NHI7vPZwxfX/KwM4/IcanOh0RqV5UJygzOwLvGtQqYCJQ/qd38J/gScC+Rg5N6qBXu+5cPvQiJqSNZ9GWr5m//lO+zVzLgh8+Y8EPn9E5qSMn9klnVN90erXrHulwRSTKRG2CMrPReMnpSeDPgRuCC8xsO16RxNbArgOpPOQnUSYhNp4T+4zgxD4jyMjdwcc/fMHCjYvYmZfFa6vf4bXV73B4h8MY1SedE/qMIKV1u0iHLCJRICoTlJn1x5t89jbn3LQqm18A7jSz5UAycDPwaCOHKPXUo21XJqSN5+Ih41iz83s+2fglX2xeyvrsTazP3sTz38zhyM5HMLL3UNJ7HkNKYvtIhywiERKVCQq4Dm/Ov/vNLPgep+nAHcDDwLd4Nxo/hdfLkiYkxhfD4C4DGNxlAFcNu5ilGSv4ZMOXLNv+LSszHSszHU8v+ScDO/+Ekb2G8tPeQ+mYmBLpsEWkEWmqI4kq+4r2s2TrCr7YspRvtq+mpOzHNaniYuJI73kMFx11Fj3bdsPn89X5fTU3oEh0OuTJYkUaS3JCG8b0G8mYfiP54PtPeWbpy5SUlQJQUlbCZ5u/4rPNX9G1TSeG9RjC8B5DGNx5AHGxNX+Uq84NWL76MaAkJRLFlKAkar266u2K5BTMh48d+3fx9ncLePu7BSTGtSat2yCO6TaYtG6D6NwmtdL+tc0NqAQlEr2UoCRq1Tx3n597T72ZJRkrWJKxgs17MvhyyzK+3LIMgO7JXRjSbSBpXQdxZJcjmtzqxyLiUYKSqFXbqsXWqT/WqT+XpJ1H5v7dLMtYyfIdq1mZ6di2L5Nt6zJ5b90n+Hw+YmNiq+2JaW5AkeimBCVRa0La+GrXlao6p1+XNqmcPmAMpw8YQ2lZKd9nbWT5jtUs376a73b/UG1yivXFcMJhx5JfXEBifOsG/7eISOiUoCRq1WfV4tiYWI7odDhHdDqci448m/ziAlbt/I631y7g20xHaWCW9VJ/GXPXvMeb7gP6pvRiUOcBDOr8EwZ26k+71m0b5d8nIrVTggqBSpUb36g+6Yd0jhPjWzM8UO0HkFu4jzW7vmf1znWs2bmu4gbh9dmbeGvtfAC6JnfmJx378JOOfRmQ2o++Kb1IiEtotj//5vrvkqZPCaqOVKrcPLRtlcyInkczoufRABSUFPLd7h9YvXMdq3d+x3e7f2DHvp3s2LeT/2z6CvCGAzskppCVn1OxztWh/PyjKSHocy3RTAmqjlSq3Dy1jmvFkK4DGdJ1IAClZaVs3rONdVkbWLf7B9ZlbWTz3oxqizWKSov4x1cvkV+cT9+U3hyW0pPWca1qPV60JQR9riWaKUHVkUqVW4bYmFj6duhF3w69GNv/RAAKigu47NXfVbt/QUkhM5e8DHj3Z3Vv24W+Kb3o26E3fVN60yelJymt21XMehFtCUGfazmYMn8Zewv3sTsvm9152WTl51R8vzs/B/Bz/cgr6dQAVbFKUHVUW8mzNG+t41vTqYaff5v4JIb3HMLG7C1s2buNjNwdZOTu4LPNS37cJyGJ3u2606td92rfAyKXEPS5brn8fj+5RfvJyd9DTsFecgr2kh34PrtgD9n5OewKJKTSaiphy/l8PnbnZStBRVJdS56learp53/V8Isrej7FpcVs3rONDTlb2JCzmQ3Zm9m0J4P9RXms2fU9a3Z9X+P7J8UnsWzbSrq06UTnNqkkxMYfNKZwXMvS57p20XS9sC7xFJUWs69wP3sL95FbtI/cwn3sLdzH3sJccvK9xJNTsJec/L3kFO6tNfEEa5vQho5JHUhN6kBqYgqpSR3oGHjs0bYrqUkdGuTfq8liQxBtH1ZpXPX5+fv9frIL9rBlzzY278lg8dblrN75HX5q/3/XoXV7urRJpXNyJ7q06UhqYkdSEtuR0tr7+jbTMXPJPw9ILJNGTGzSRRvRpOr1Qqj/OQ6F3++noKSQ/cV55BXle4/FBSzNWMGH6z+j1P9jUvHho3NSR/w+r0K1oKQwpGO1iU8kJbE9Ka3b0aG191j+vGNiCp0CiSghLiHc/8wKtU0WqwQl0sjKE8KuvCzat2pLWrdBtIpNYMf+nWTu282uvKyK+7VClRjXmnMHnkZyQhvaJCSRnNCG5ISkwPdJJMYnEhcTG+Z/UfhEU7K89s3bqh3+7JTUkSfGTa7UVlpWSmFJEQUlhRSUFAQeD/JVHEhCxfkViWh/cT55xfnU9/dybEws7RKSadsqmbat2tC2VXLF8w6J7UipkoTq0lOvSbh+VprNXCSKHOzertKyUrLyc8jcv5vMfbvI3L+b7PyciusEOQV7ycrPqfa1+SUF/HPlm7UePy4mjsS4VrSOa0Xr+Na0jmtFYlzrwPNWlZ63imtFfEwc8bHxgcc44mLiSIiNJy7wPD4mjrjYOOJ8sfh8PmJ8McT4Yip9H0P59z58QY/4/ZThB7+fTzct5pml/6SotBgor3B8kfziAkb2GlqxXxl+ysrKKPGXUloW+PKXUVpWSklZKWV+77HUX0ppWVngsXxbmbctsL24tITismLvsbSY4rIfH2u6XrgrL4vr5t3u7VdaTFHgNeHUOq4VSfGJtIlPJCkhiTbxiSzdtrLG/aedfQ9tWyWTGNc6pGVo6quxqlGVoESiTGxMLJ3bpNK5TSpHdjmi2n2uffPP7MrLPqC9TXwip/1kNPuK8thflMe+ov3sD/o+v6SQkrIScotKyC3a39D/lENWVFrMzCWzmblkdqRDqWTn/t2VnvvweQk++CuQ7FtVba/4ak2bhMQDElFSfCKx1fRya+vRdU3u3GD/1uo0VjWqEpRIEzQh7bwaijZ+UesvCL/fT3FZCQXF3jBUfmA4Kr/4wKGp/OICCksKvV5FWQklpSUUlRVTUlpSqadR3l5WVub1bvxllPnL8Pv9lR7Lyns/Qdt9Ph8+fPh8vlp7IW1bJeMDfL4YfECsL5bYmBhiY2KJ88USE3iMjQm0V3wfS6wveL+YH/fzxXo9wIreYeXH77M28sH6TystmhkfE8fPjxrHyN5DiY+NJyEm3nuMjW/wnks0FbQ01u0JSlAiTVB95ikEryQ4IfALtR3RNedgKNd8GsOYfiMZ0Klf1FwTq+/PvCE01u0JSlAiTdShzlMYbaKph1Au2s5xtMTTWD8rJSgRiQrR1EOQ2jXWz0oJSkSiRrT0EOTgGuNnFdOg7y4iIlJPSlAiIhKVlKBERCQqKUGJiEhUaolFErEA27dvj3QcIiItXtDv4gOmz2iJCao7wMSJEyMdh4iI/Kg7UGlNmpaYoBYDo4BtQN0WQxERkYYSi5ecFlfd0OKW2xARkaZBRRIiIhKVlKBERCQqKUGJiEhUUoISEZGopAQlIiJRSQlKRESikhKUiIhEJSUoERGJSi1xJol6M7OjgSeBNGA9cJVz7oC7nyPJzE4DHgAGAJnAQ865GWbWCsgFioJ2/8w5918RCLOCmV0FzAAKg5qvA2YDjwMX4c348Yhz7v7Gj7AyM5uIF2+wRGA+MI4oO8dmlg7Mc851CTxPoJbzamY/B+7Du7P/Y+AK51xmBOPtAjwKnAr4gLeBG51z2YHtzwM/B0qC3ibNObc+gjHX+n8tCs/xviq7xAGtgJ7OuYxInmMlqDoK/MeeC0wFRgMXAu+ZWR/n3N6IBhdgZr2BOcDleLEOB941sw3AbiDLOdctYgFWbxjwsHPuT8GNZnY/YEB/oD3wjpltdc49H4EYKzjnZgGzyp+b2VDgPeAPwBCi5BybmQ/4b2BKlU13U8N5NbPBwNPAmcBXwIPAy8ApEYx3JrAH6AfEAy8A04FLAtuHAec5595p6BirqiXmGj8H0XiOnXPJQfvEAQuAj5xzGYHmiJ1jDfHV3UlAvHNuqnOu2Dn3MvAtcHFkw6qkL/CSc+4151xZoHf3EXACXrL6OoKx1aSmuC4HJjvnsp1zG/D+U01qzMAOxszi8ZLVXc65b4iuc3w3cA3w1yrttZ3XXwJvOuc+dc4VALcCJ5jZgEjEa2YxQBlwt3Nuv3MuB/gHcGJgeyIwkMid85rOcW2fg6g6x9W4Be8PgTsh8udYCaruBgOrq7StwftrKSo45xY6535d/tzMOuJNjLsM76+gLma23Mx2mNn/M7OekYo1EF8s3nDppWaWYWbrzOxPZtYBb/hjVdDuUXWuA64D8oEnAs+j6Rw/6ZwbjvdXOgBmlkLt53Vw8DbnXB6wmcY57wfEG/gj6zzn3Lqg/c7D+zwDHIM37PQPM9tpZkvN7JxGiLXcATEH1PY5iKpzHMzMegB/Bn7tnCsLNEf0HCtB1V0ykFelLQ9IikAsB2Vm7YE3gC/xhvv2A//BG8s3vF+sr0UsQE9nvP8s/4c3hHMR3l941we2B5/vqDrXgSHfP+D1nspnXI6acxw0PBOsfCinpvMasc94DfFWYmY34yWoWwJNbYGFeD2DHsBk4F+Ba8UNrpaYa/scRPM5/h3wjnMuuLcU0XOsa1B1tx/vYniwJKDqBcaIM7Mj8JLSKmBi4K+hm6rscxOw08x6O+c2RyBMnHPbgTFBTV+b2TS88XmofL6j7VyfgTf89FZ5g3Mu6s5xFfsDjzWd16j8jAeGUqfhFaGc4pxbA+Ccew/v+l+5OWZ2JXAu8E2jBxpQ2+eA6D3HsXjDv5UWyov0OVYPqu5W4f01FGwglYdLIs7MRuP1ml4HLgqMc2Nm95jZoKBdEwKPBY0cYgUzO9LM7q7SnIAX03Yqn+9oO9fjgX8FDYVE5TkOFqh8q+28VvqMm1kScBgRPO9m1hZ4HxgBpAf/dW9m48zs8iovKf/8RMxBPgdRd44Djg88zg9ujPQ5Vg+q7hYAPjP7HV6Z7oV4108iPUxWwcz6A/OA25xz06psTgOONbPy6qdHgbecczsbM8YqcoDfm9kWvMqmocANwG/wClDuNLPleMMiN+PFHC1GAn+p0haN57iqF6j5vL4EfGpmJwGfA/cDy5xzayMRaMDLeH9IjwpcrwkWCzxqZquBJXgFS8cDVzduiAeo8XNgZtF4jsH7PH8R/AdXQETPsXpQdeScK8IberoQyAJuwyu9jKZfPtfhjRnfb2b7gr4exCsvzQbWARvw7tG4NGKRAs65rXhDBZOAvXgl8vc6514B7gBW4iWqxYFtT0Yo1Or0BaqO6UfdOa5GjefVObcCuCrwfBdwJPCzyIQJZpYGnAWkA5lBn+ctgXhfx/t/OBvv8/N74Bzn3KZIxRxQ4+cg2s5xkL4c+HmO+DnWiroiIhKV1IMSEZGopAQlIiJRSQlKRESikhKUiIhEJSUoERGJSkpQIiISlXSjrkiIzOw5vGlhanI33izyC4C2zrlGmcYmMF3Nf4DL6nPjp5n5gXHOuXl12PdxYLFz7v9Cj1SkbtSDEgndjXizgnfHW4YFvJtJy9umAJ8Fvt9fzesbyg3AN4cwK0F3vGmF6uIe4B4zS63nsUQOSjfqihwCMzsKWAH0C6yvFKk4WgOb8CZTXdlIx3wG2OScu6sxjictj4b4RBpAYK61iiG+wPDZBLwF6gxvmZFf4i3ZcSneNDK3OudeCLy+LfAw3hIkfuBDvKXOa1oy4RdATnlyMrO+wA94U0k9AvQCPsBbzmQK3szg24BrAzNWVxriM7OP8JYjPxr4L7w1ix5yzs0MOuarwNNmNtk5V1zvkyVSAw3xiTSeB4Df4k3MeRiwFC8xjcD7ZT/DzMrXbHoKL5GdjrckiR94N7Akd3XOBqpbkvtevOXRT8MbjlyON/w4HG/hv5nVvKbcLXhDfkPxktvfzSx4GfMPgNTAe4mEnRKUSOOZ7pxbEFgyYh7eGkB/ds45vF5OItDPzA7H6xFd4pxbHOgVXYo3oecZNbz3sXgTwFY1OfAen+ItPLfKOfdYYE2l6UDvQG+tOh8556YH4rsVb8QlrXxjYCmX9YFji4SdEpRI4wleujwP2BC0Gm/5+jqt8JYFB3DlM3gDu4E2HLgmWbmueLNjH+yY64OeBx+zOhXFFs65vYFv46vssxvoUsPrRQ6JrkGJNJ6q12mqrr1TLi6w71C8ob1gWTW8pgzwHcIxq1NUTVvVY8QCpSG8p0idqQclEn1W4/VU2jjn1jnn1uEVNDwEHFHDa7YDnRspvmCdAscWCTslKJEoE7jm8wbwvJmNMrOBwPN4xRVranjZEryKu0ZjZu2BPngLH4qEnRKUSHS6HK8U/XW8BNAeOM05l1PD/m/hVfs1phPxek/LGvm40kLoRl2RZsDMkvCWFz/DObe0kY45G68q8N7GOJ60POpBiTQDzrk8vGtU1zXG8cysO16PbXpjHE9aJiUokebjb0CamdVUih5OtwO3O+dqqioUOWQa4hMRkaikHpSIiEQlJSgREYlKSlAiIhKVlKBERCQqKUGJiEhU+v/RtT2C8cH46AAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Solution\n",
    "\n",
    "plot(results.I, 'g-', label='simulation')\n",
    "plot(data.insulin, 'go', label='insulin data')\n",
    "\n",
    "decorate(xlabel='Time (min)',\n",
    "         ylabel='Concentration ($\\mu$U/mL)')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Exercise:** Using the best parameters, estimate the sensitivity to glucose of the first and second phase pancreatic responsivity:\n",
    "\n",
    "$ \\phi_1 = \\frac{I_{max} - I_b}{k (G_0 - G_b)} $\n",
    "\n",
    "$ \\phi_2 = \\gamma \\times 10^4 $\n",
    "\n",
    "For $G_0$, use the best estimate from the glucose model, 290.  For $G_b$ and $I_b$, use the inital measurements from the data.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Solution\n",
    "\n",
    "I0, k, gamma, G_T = best_params"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(130, 11)"
      ]
     },
     "execution_count": 18,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Solution\n",
    "\n",
    "I_max = data.insulin.max()\n",
    "Ib = data.insulin[0]\n",
    "I_max, Ib"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(289, 92)"
      ]
     },
     "execution_count": 19,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Solution\n",
    "\n",
    "# The value of G0 is the best estimate from the glucose model\n",
    "G0 = 289\n",
    "Gb = data.glucose[0]\n",
    "G0, Gb"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "18.455830724497748"
      ]
     },
     "execution_count": 20,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Solution\n",
    "\n",
    "phi_1 = (I_max - Ib) / k / (G0 - Gb)\n",
    "phi_1"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.8184889071339871"
      ]
     },
     "execution_count": 21,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Solution\n",
    "\n",
    "phi_2 = gamma * 1e4\n",
    "phi_2"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}
